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ABSTRACT 

We  state  and  prove  a  useful  theorem  on  manipulating  rotation  order  which, 
while  not  new,  is  barely  present  in  the  literature.  This  theorem  allows  the 
order  of  a  sequence  of  rotations  to  be  reversed,  provided  that  the  sense  of  the 
axes  of  rotation  is  changed  from  “body”  to  “space  fixed”  or  vice  versa.  We  use 
the  theorem  to  aid  calculations  in  geodesy  (constructing  a  local  north-east- 
down  coordinate  system)  and  aerospace  theory  (relating  yaw-pitch-roll  rates 
to  vehicle  angular  velocity) .  The  new  notation  here  sheds  light  generally  on  the 
held  of  orientation  theory,  as  well  as  giving  insight  to  standard  terms  relating 
to  wind  direction  used  for  treating  ship  motion.  Although  we  present  our 
analyses  in  the  style  of  a  tutorial  in  the  general  subject  of  spatial  orientation 
theory,  there  is  new  notation  here,  along  with  alternative  and  novel  ways  of 
treating  problems  that  are  often  seen  as  difficult  or  obscure  by  practitioners. 
This  report  follows  on  from  the  2005  DSTO  report  DSTO-TN-0640,  but  is 
completely  self  contained,  and  DSTO-TN-0640  need  not  be  read  beforehand. 
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A  Pseudo-Reversing  Theorem  for  Rotation  and  its 
Application  to  Orientation  Theory 

Executive  Summary 

This  report  is  a  much-evolved  follow-on  from  the  2005  DSTO  publication  DSTO-TN-0640, 
“Using  Rotations  to  Build  Aerospace  Coordinate  Systems”,  that  explained  the  construction 
of  coordinate  systems  used  in  aerospace  calculations.  That  report  followed  a  step-by- 
step  approach  to  implement  its  calculations.  In  the  current  report  we  rephrase  those 
calculations  in  a  more  efficient  language,  while  incorporating  a  very  useful  theorem  that 
is  known  but  almost  absent  from  the  literature.  This  “Pseudo- Reversing  Theorem”  allows 
the  order  of  a  sequence  of  rotations  to  be  reversed,  provided  that  the  sense  of  the  axes  of 
rotation  is  changed  from  “body”  to  “space  fixed”  or  vice  versa.  The  current  report  is  self 
contained,  so  that  familiarity  with  the  content  of  DSTO-TN-0640  is  not  necessary. 

The  current  report  places  the  theorem  and  the  reworked  examples  of  DSTO-TN-0640 
into  the  greater  context  of  orientation/rotation  theory.  We  first  introduce  the  theorem, 
then  establish  a  solid  mathematical  language  necessary  for  quantifying  the  orientation  of 
an  object.  We  cover  the  background  of  how  to  rotate  a  vector,  using  either  a  matrix 
or  a  quaternion.  We  then  rework  the  examples  in  DSTO-TN-0640:  constructing  a  local 
north-east-down  set  of  axes  from  a  given  latitude  and  longitude,  and  calculating  where  a 
pilot  must  look  to  see  a  distant  aircraft.  We  also  make  an  extended  revisit  to  the  subject 
of  conversions  within  the  Distributed  Interactive  Simulation  environment  for  handling 
orientation  information,  since  this  often  causes  problems  to  practitioners  who  must  deal 
with  several  coordinate  systems  at  once.  We  end  the  main  report  by  showing  how  the 
Pseudo-Reversing  Theorem  can  be  used  to  simplify  some  of  the  concepts  behind  dead 
reckoning  an  object’s  changing  orientation. 

The  report  ends  with  an  appendix  that  applies  its  notation  and  general  approach  to 
the  task  of  constructing  the  appropriate  course  a  ship  must  steer  in  order  for  the  wind 
to  appear  to  come  from  some  given  direction  with  some  given  speed.  This  is  a  nontrivial 
problem  that  is  handled  well  in  a  novel  way  by  the  orientation-matrix  language  of  this 
report,  although  its  solution  doesn’t  require  the  Pseudo-Reversing  Theorem. 
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1  Introduction 

Since  publishing  a  DSTO  report  [1]  covering  the  basics  of  rotation/orientation  theory 
in  2005,  I  have  received  several  emails  from  practitioners  with  various  questions  prompted 
by  its  subject  matter.  These  questions  suggest  to  me  that  a  second  edition  of  that  report 
would  be  useful. 

There  are  several  reasons  to  update  the  2005  report.  The  subject  of  6  degree-of-freedom 
modelling  is  a  comparatively  modern  one,  requiring  the  use  of  fast  computers  to  update  an 
object’s  state  sufficiently  accurately  to  represent  the  real  world  usefully.  Perhaps  because 
of  this,  straightforward  information  about  the  subject  has  yet  to  appear  in  nonspecialised 
textbooks.  The  calculations  in  [1]  were  written  to  fill  a  perceived  gap,  having  both  ped¬ 
agogical  and  logical  content.  But  they  were  certainly  longer  than  necessary.  Only  after 
writing  that  report  did  I  become  more  aware  of  a  known  theorem  that  would  have  simpli¬ 
fied  some  of  its  analysis.  I  think  that  this  theorem  is  not  given  the  significance  it  deserves 
in  the  literature;  in  fact,  while  some  practitioners  are  aware  of  it,  a  statement  of  it  with 
a  discussion  of  its  use  is  hard  to  find  anywhere — perhaps  because  having  no  widely  used 
name  makes  it  difficult  to  search  for.  I  have  called  it  the  Pseudo- Reversing  Theorem  here, 
or  PRT.  It  has  been  called  the  Rodrigues  Transposition  Theorem  in  [2],  but  no  provenance 
is  given  there  for  the  choice  of  name,  which  doesn’t  appear  to  be  used  anywhere  else.  An 
example  of  the  theorem  is  given  in  [3],  which  gives  a  proof  in  the  context  of  orthogonal 
axes,  although  the  theorem  doesn’t  actually  require  such  axes. 

The  Pseudo- Reversing  Theorem  can  often  be  invoked  to  give  a  different  pedagogical 
basis  to  the  many  analyses  and  recipes  abounding  in  orientation  theory  that  can,  for  some, 
seem  quite  opaque.  Part  of  the  subject’s  difficulty  is  due  to  a  divide  in  the  orientation 
community:  some  books  base  their  analysis  on  the  concept  of  an  “active  rotation”,  whereas 
others  use  a  “passive  rotation”.  I  suggest  that  knowledge  of  the  Pseudo- Reversing  Theorem 
forms  a  good  bridge  between  these  two  approaches,  so  that  one  can  more  easily  appreciate 
the  reason  why  both  approaches  have  historically  been  used. 

Confusion  over  whether  a  sequence  of  rotations  is  active  or  passive  can  result  in  that 
sequence  being  written  in  the  wrong  order,  possibly  with  wrong  signs  of  the  angles  turned 
through.  But  orientation  theory  is  a  subject  in  which  a  procedure’s  lack  of  validity  will 
quickly  become  obvious  when  it’s  implemented  on  a  computer.  It  does  not  seem  to  me  that 
a  great  deal  of  pedagogical  effort  has  found  its  way  into  the  proofs  of  some  of  the  subject’s 
recipes;  and  when  they  work  well,  there  is  little  incentive  for  anyone  to  highlight  the  places 
where  the  analyses  are  unclear  and  to  provide  an  explanation  for  why,  nonetheless,  those 
analyses  work.  I  have  tried  to  address  that  in  this  report  by  explaining  some  representative 
analyses  in  detail.  And  although  I  use  only  active  rotations  here  and  in  [1],  I  do  discuss  how 
passive  rotations  can  be  used,  and  where  they  fit  into  the  general  scheme  of  the  subject. 

Another  reason  for  this  follow-up  report  is  to  establish  a  good  set  of  notation  that  helps 
simplify  rotation/orientation  analysis.  I  rework  the  examples  in  [1]  using  this  notation. 
I  have  added  an  appendix  that  uses  the  notation  to  aid  calculations  of  wind  velocity 
as  perceived  by  a  ship.  This  is  a  standard  task  in  navigation  that  is  rendered  more 
transparent  by  an  appreciation  of  the  difference  between  proper  vectors  and  coordinate 
vectors,  as  discussed  in  Section  3. 

Finally,  Section  5.4  on  “DIS  conversion”  was  prompted  by  a  question  that  I’ve  been 
asked  several  times  regarding  the  DIS  example  in  Section  4.3  of  [1],  an  example  that  has 
been  useful  for  DIS  practitioners.  That  original  calculation  was  certainly  correct,  but 
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readers  were  apt  to  confuse  the  roles  of  Earth-centred  Earth-fixed  and  local  north-east- 
down  and  so  arrive  at  a  wrong  answer.  To  address  this,  I’ve  rewritten  the  discussion  of  [1] 
differently  and  in  far  more  detail  here. 


2  The  Pseudo-Reversing  Theorem 

The  following  thought  experiment  encapsulates  the  core  of  the  calculations  in  this  report. 

Consider  the  viewpoint  of  an  audience  watching  the  famous  mime  Marcel  Marceau  on 
stage.  Marcel  orientates  his  body  in  some  way,  and  then  rotates  his  hand  about  his  wrist. 
We,  the  audience,  seek  to  describe  that  hand’s  final  spatial  orientation  from  a  knowledge 
of  these  two  procedures.  Orientating  Marcel’s  body  can  be  described  by  some  action,  an 
“operator”,  which  we  write  as  O  and  assume  known.  To  incorporate  the  new  orientation 
of  Marcel’s  rotated  hand,  we  can  write  another  operator  that  rotates  that  hand  about  the 
wrist.  But  simply  saying  “the  rotation  was  about  the  wrist”  is  useful  mathematically  only 
if  we  know  where  the  wrist  currently  points.  So  we  wish  to  specify  the  vector,  independent 
of  Marcel’s  body,  along  which  his  wrist  points.  The  act  of  bodily  orientating  himself  has 
changed  the  direction  of  Marcel’s  wrist,  in  which  case  we  must  distinguish  between  the 
initial  direction  (vector)  of  his  wrist,  here  denoted  by  the  subscript  “wrist”,  and  the  final 
direction  vector  of  his  wrist,  here  denoted  by  the  subscript  “(wrist)”. 

If  we  write  the  operation  of  Marcel  orientating  himself  bodily  as  O,  and  follow  this 
with  a  rotation  -R{wriat)  of  his  hand  about  the  latest  direction  of  his  wrist,  then  the  entire 
procedure  of  orientating  Marcel’s  hand  is  written 

final  orientation  =  O  — >  R(  wrist)  ■  (2-1) 

Now  realise  that  Marcel  can  achieve  exactly  the  same  result  by  first  rotating  his  hand 
about  his  wrist,  and  then  orientating  himself  bodily  as  before.  We  still  describe  the  bulk 
orientation  with  O,  but  we  must  now  rotate  the  hand  about  a  different  vector:  the  initial 
direction  of  Marcel’s  wrist.  Write  this  rotation  as  RWI ist  to  obtain 

same  final  orientation  =  Rwris t  — ►  O  .  (2-2) 

We  have  “almost”  swapped  the  two  procedures  in  the  mathematical  description  of  Marcel’s 
movements,  provided  we  understand  the  operator  that  rotates  the  hand  to  be  different  in 
each  description.  This  simple  result  forms  the  core  of  what  we’ll  soon  call  the  Pseudo- 
Reversing  Theorem: 

O  >  -^(wrist)  -^wrist  *  O  •  (2-3) 

In  the  descriptions  to  follow,  we’ll  use  a  phrase  such  as  “the  latest  snapshot  of  the  vector 
describing  Marcel’s  wrist”  to  reinforce  the  fact  that  there  can  be  two  wrist  vectors  being 
discussed,  which  refer  to  the  direction  of  Marcel’s  wrist  at  different  times  in  an  orientation 
procedure.  The  vector  “wrist”  describes  where  his  wrist  was  initially  and  is  unchanging;  the 
vector  “(wrist)”  describes  where  his  wrist  is  now,  and  is  subject  to  change  as  he  orientates 
his  body.  We  might  consider  the  initial  vector  as  being  rotated  to  become  the  final  vector, 
but  it’s  useful  to  view  the  initial  vector  as  set  in  stone,  and  a  snapshot  of  this  vector  is 
then  physically  rotated  to  become  the  final  vector. 

More  generally,  consider  a  sequence  S  of  rotations  about  the  latest  snapshots  of  vec¬ 
tors  Ui,u2,u3.  We  wish  to  describe  the  following  set  of  rotations.  First,  rotate  snapshots 
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of  u3,u2,u3  around  u3  to  give  a  new  set  labelled  ufiu'2,u'3  respectively.  Now  rotate 
snapshots  of  each  of  these  vectors  around  u2  to  give  a  new  set  labelled  u" .  u2.  u3  respec¬ 
tively.  Now  rotate  snapshots  of  each  of  these  vectors  around  u3  to  give  a  new  set  labelled 
u'C,  u'2  ,  u3  respectively.  Finally,  rotate  snapshots  each  of  these  vectors  around  u'fi  to  give 
a  new  set  labelled  Vi,v2,v3  respectively.  We  can  write  the  sequence  of  rotations  between 
initial  and  final  vectors  as 

{wl5  u2,  u3}  ->  RUi  -*•  Rui2  -*•  Run  ->  Rum  ->  {Wi,  V2,  U3}  .  (2.4) 

S  is  the  sequence  of  four  rotations  taking  initial  to  final  vectors.  But  this  linguistic  and 
notational  description  of  the  sequence  is  tedious — and  we  almost  certainly  will  not  need 
the  intermediate  vectors — so  instead  we’ll  describe  the  same  procedure  as  follows.  Rotate 
snapshots  of  u3,u2,  u3  first  around  u3,  then  around  the  latest  snapshot  of  u2,  then  around 
the  latest  snapshot  of  u3,  and  finally  around  the  latest  snapshot  of  Ui,  to  give  a  new  set 
Vi,v2,v3.  We  will  write  this  as 

{u1,u2,u3}  — >  R±  — >  R( 2)  — >  R( 3)  — >■  R{  1)  — >  {v1,v2,v3}  .  (2.5) 

The  above  discussion  of  Marcel  Marceau  encapsulated  in  (2.3)  allows  us  to  replace  the 
Ri  — >  R( 2)  in  (2.5)  with  R2  —>  R1.  This  doesn’t  affect  the  sense  of  the  next  vectors  (3) 
and  (1),  since  the  procedures  R3  —>  R (2)  and  R2  —■ ►  R\  have  the  same  effect — they  produce 
the  same  final  orientation.  In  that  case,  (2.5)  becomes 

S  =  R-2  — >  R\  —>  R( 3}  — ►  R{ i)  • 

Now  again  use  the  same  logic  of  applying  (2.3):  think  of  R2 
and  swap  it  with  R,/3\ ,  remembering  to  change  R(3\  to  R3: 

S  =  R3  R2  — >  Ri  — >  R(i)  ■ 

Finally  apply  the  same  idea  again  to  arrive  at 

S  =  R3  -*•  R3  ->  R2  ->  R3  .  (2.8) 

Notice  that  we  cannot  combine  in  (2.7)  into  one  rotation:  R1  refers  to  a  rotation 

around  ult  while  Rn\  refers  to  a  rotation  around  the  very  latest  “incarnation”  of  ult  which  in 
general  is  a  completely  different  vector  to 

An  inductive  argument  should  make  it  clear  that  the  same  reversing  procedure  holds 
regardless  of  the  length  of  the  sequence.  There  can  be  repeated  subscripts,  and  the  initial 
set  of  vectors  need  not  be  mutually  orthogonal.  For  any  number  of  rotations,  then,  we 
have  proved  the  Pseudo -Reversing  Theorem: 

R3  — >  R( 2)  — ►  R( 3)  — ►  R( 4)  —>•••■  =  •  ■  ■  — >  ^4  — R3  ~ *  R2  y  Rl  ■  (2-9) 

We  can  equally  well  express  the  right-hand  side  of  (2.9)  as  a  purely  operator  expression, 
in  which  case  it  uses  no  arrows  because  operators  act  from  right  to  left: 

(2.10) 

Equation  (2.10)  is  a  very  useful  statement  of  the  Pseudo-Reversing  Theorem.  It  says 
that  when  rotating  about  latest  axes,  we  need  only  write  down — from  left  to  right — the 
sequence  of  rotations  nominally  around  latest  axes,  but  writing  them  as  operators  that 
rotate  about  space- fixed  axes,  which  will  then  be  automatically  understood  to  act  from 
right  to  left.  We’ll  do  this  in  the  examples  to  follow. 


— >  R 


<2> 


R 


<3> 


R 


<4> 


—  Ri  R2  ^3  -R4 


(2.6) 

R3  as  the  orientation  O, 
(2.7) 
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2.1  A  Related  Scenario  in  Aerospace  Kinematics 

The  above  discussion  of  Marcel  Marceau  is  actually  central  to  quantifying  aerospace  kine¬ 
matics,  and  in  the  interests  of  clarity,  we  will  reword  it  using  the  more  complex  language 
of  aircraft  axes. 

Suppose  an  aircraft  is  flying  straight  and  level  with  its  nose  and  wings  pointing  in 
some  given  directions  relative  to  north,  east,  and  down.  The  pilot  executes  three  quick 
manoeuvres.  First,  a  yaw  through  10°,  followed  by  a  pitch  around  the  wing  through  20°, 
followed  by  a  roll  around  the  nose  through  30°.  Assuming  the  directions  of  north,  east, 
and  down  haven’t  changed  appreciably  during  the  manoeuvres,  what  is  the  final  aircraft 
orientation  relative  to  north,  east,  and  down? 

It  suffices  to  specify  this  orientation  by  giving  coordinates  that  quantify  the  aircraft’s 
three  “body”  basis  vectors  for  its  nose,  starboard  wing,  and  onboard-down  direction.  These 
vectors  can  be  envisaged  as  arrows  embedded  in  the  aircraft’s  body.  Call  these  pre¬ 
manoeuvre  vectors  “nose”,  “starboard  wing”,  and  “down”,  n0,s0,d0,  respectively.  Their 
coordinates  relative  to  north,  east,  and  down  are  specified  in  the  scenario.  We  will  be  extra 
clear  here  by  using  the  above  laboured  description  of  intermediate  rotations  once  more. 
The  first  manoeuvre  rotates  these  vectors  around  d0  to  produce  n1,s1,d1.  (Of  course, 
d1  =  o?o>  but  a  good  naming  convention  helps  prevent  confusion.)  The  second  manoeuvre 
rotates  n1,  s1,  dx  around  S}  to  produce  n2,s2,d2.  The  third  manoeuvre  rotates  n2,s2,  d2 
around  n2  to  produce  n3,s3,d3.  We  require  the  coordinates  of  n3,  s3,  d3  relative  to  north, 
east,  and  down.1  This  sequence  of  three  rotations  can  be  written 

{n0,  s0,  d0}  —>■  Rdo(W°)  — y  RSl(20°)  ^n2(30°)  — *  {n3,  s3,  d3}  .  (2.11) 

Each  of  these  rotations  requires  computational  effort  because  each  rotation  (except  the 
first)  requires  a  vector  produced  by  the  previous  rotation,  and  so  cannot  be  constructed 
in  advance  of  the  manoeuvres;  the  intermediate  vectors  must  be  calculated.  This  sort  of 
calculation  was  done  in  [1],  Now,  however,  we  can  use  the  PRT  to  write  the  final  result 
as  the  simpler  sequence  of  three  rotations  solely  around  the  initial  basis  vectors: 

{n0,s0,d0}  -  Rno( 30°)  -  RSo( 20°)  -  i^(10°)  -  {n3,s3,d3}.  (2.12) 

This  sequence  is  easier  to  implement  than  (2.11)  because  no  intermediate  basis  vectors 
need  be  calculated. 

The  Pseudo-Reversing  Theorem  has  connected  these  two  descriptions  of  a  sequence  of 
rotations  of  the  aircraft.  Both  rotation  sequences  are  necessary  and  useful  for  application 
in  geodesy  and  flight  modelling.  Rotations  about  carried-along  axes  closely  match  our 
physical  experience:  they  are  what  a  pilot  flies;  but  such  rotations  are  not  always  econom¬ 
ical  to  describe  mathematically  because  of  the  need  to  calculate  intermediate  vectors  to 
implement  the  rotations.  Rotations  about  space- fixed  axes  are  mathematically  economical 
because  intermediate  vectors  need  not  be  calculated,  but  these  rotations  do  not  match  our 
physical  experience,  which  can  make  them  more  difficult  to  visualise.  (For  example,  a  pilot 
performing  aerial  manoeuvres  doesn’t  concentrate  on  rotating  the  aircraft  about  the  east 

1In  practice,  the  final  directions  of  north,  east,  and  down  can  change  during  the  scenario,  but  this  has 
no  real  bearing  on  the  discussion;  we  need  only  choose  a  coordinate  system  to  express  everything  in,  and 
the  original  north-east-down  directions  suffice  as  its  axes.  We  could  easily  accommodate  any  change  in 
the  final  directions  of  north,  east,  and  down  using  Section  3. 
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direction;  that  would  require  a  very  complex  interaction  with  the  flight  controls.)  The 
PRT  allows  the  best  of  both  worlds:  easy  mathematics  and  easy  visualisation.  It  says  that 
any  sequence  of  rotations  in  one  description  is  converted  to  the  equivalent  sequence  in  the 
other  by  reversing  the  order  of  the  rotations  and  changing  their  sense  to  that  of  the  other 
description. 


3  Language  of  Vectors  and  Rotations 

Like  its  position,  the  orientation  of  an  object  must  be  specified  relative  to  some  unchanging 
reference  orientation.  Euler’s  theorem  of  rotation  states  that  any  orientation  of  an  object 
can  always  be  reached  by  applying  a  single  rotation  to  that  base  orientation  [4].  Hence 
we  can  quantify  any  orientation  by  writing  down  this  single  rotation.  Alternatively,  we 
might  prefer  to  restrict  all  rotations  to  a  prescribed  set  of  axes,  in  which  case  more  than 
one  rotation  will  generally  be  needed,  and  we’ll  have  to  specify  the  order  in  which  the 
rotations  are  performed. 

These  procedures  necessitate  rotating  one  vector  about  another,  together  with  good 
notation  that  keeps  track  of  the  various  coordinate  sets  that  might  be  involved.  We  will 
begin  with  a  proper  vector  v,  being  an  arrow  that  signifies  a  directed  number  such  as 
velocity.  Note  the  distinction  between  a  proper  vector — an  arrow  that  might  be  attached 
to  one  point  in  space,  such  as  those  comprising  the  velocity  field  of  a  moving  fluid,  and  a 
position  vector  (or  radius  vector ),  which  is  an  arrow  whose  tail  is  anchored  to  the  origin, 
with  its  head  serving  to  quantify  a  point  in  space.  For  example,  one  position  vector  minus 
another  position  vector  gives  a  proper  vector.  We’ll  consider  this  distinction  further  in 
Section  3.2.4. 

A  coordinate  vector  [u]g  is  a  column  of  numbers  that  quantifies  v  as  a  linear  com¬ 
bination  of  a  set  of  basis  (proper)  vectors  S.  A  proper  vector  can  be  described  using 
any  one  of  an  infinite  number  of  bases  S,  S', ,  being  represented  by  coordinate  vectors 
[u]5,  [u]<y/, . . .  respectively.  These  coordinate  vectors  are  all  distinct  columns  of  numbers, 
each  associated  with  its  basis,  but  all  describe  the  same  proper  vector  v.  Proper  vectors 
can  be  manipulated  in  the  usual  way  by  “adding”  arrows,  but  a  numerical  representation 
demands  a  basis,  with  each  proper  vector  represented  by  its  appropriate  coordinate  vector 
in  that  basis.  To  avoid  notational  clutter,  we’ll  use  e.g.  S  to  refer  to  a  set  of  axes,  and  to 
the  basis  vectors  of  that  set. 

3.1  Notation  of  Coordinate  Vectors 

Orientation  analyses  simplify  when  the  basis  vectors  are  orthonormal  (mutually  orthog¬ 
onal  and  of  unit  length).  Although  the  PRT  applies  also  to  non-orthogonal  axes,  our 
geodesy /aerospace  focus  assumes  orthonormal  bases  S  and  S',  with  basis  vectors  writ¬ 
ten  as 

S  basis  =  { ex ,  ey ,  ez }  ,  S'  basis  ={  ex, ,  ey, ,  ez, }  ,  (3-1) 

and  so  on.  Now,  because  a  proper  vector  can  be  written  as  a  linear  combination  of  basis 
vectors 


V  VXGX  T  VySy  ~f“  VZSZ 

=  vxtexi  +  vyieyi  +  vziezi ,  (3.2) 
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the  coordinate  vectors  with  respect  to  bases  S,  S'  are 
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Central  to  orientation  theory  is  the  description  of  one  coordinate  system’s  orientation 
relative  to  another.  For  example,  we  might  require  the  orientation  of  an  aircraft  in  the 
Earth-centred  Earth-fixed  coordinate  system  (ECEF)  defined  in  Section  5.2.  The  position 
of  a  point  such  as  the  aircraft  axes’  origin  in  these  ECEF  coordinates  is  not  relevant  to 
specifying  its  orientation,  so  the  relevant  proper  vectors  can  be  envisaged  as  attached  to 
the  aircraft,  and  are  unchanged  if  the  aircraft  is  merely  translated.  The  orientation  of  S 
in  S'  is  quantified  by  an  orientation  matrix  n whose  columns  are  the  S'  coordinate 
vectors  of  the  S  basis  vectors:2 
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Because  basis  vectors  suffice  to  quantify  any  other  vector,  the  matrix  suffices  to  quan¬ 
tify  the  orientation  of  an  object.  The  dot  products  in  (3.4)  equal  the  cosines  of  the  angles 
between  the  various  basis  vectors,  so  /i|,  is  also  called  the  direction  cosine  matrix  relat¬ 
ing  S  and  S'.  Not  surprisingly,  /i|  is  the  identity  matrix.  Also,  some  straightforward 
linear  algebra  can  be  used  to  show  that 

tfs  =  (a4')*  =  040  1>  (3 *-5) 

where  the  superscript  t  denotes  the  matrix  transpose.  When  the  transpose  of  a  matrix 
equals  its  inverse,  it’s  known  as  orthogonal ;  its  rows  will  be  mutually  orthonormal,  and  so 
will  its  columns. 

The  mathematics  of  orientation  theory  can  be  constructed  from  the  answers  to  two 
fundamental  questions. 


First  Fundamental  Question:  Given  [v]s,  what  is  [v] ? 
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(3.6) 


So  the  orientation  matrix  transforms  one  coordinate  vector  to  another.  Equation  (3.6)  is 
an  important  equation  of  this  report,  and  we  highlight  it  by  putting  it  in  a  box: 


(3.7) 


2The  symbol  “=”  denotes  a  definition:  a  =  b  means  that  a  is  defined  to  equal  b.  We  will  occasionally 

use  the  “=”  symbol  last  in  a  string  of  equalities,  in  which  case  we  mean  that  the  last  quantity  is  being 

defined. 


Ns'  =  4  Ms  ■ 
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Figure  1:  Rotating  a  vector  u  in  a  right-handed  sense  around  the  vector  n 
to  produce  v. 

Notice  that 

[«U  =  hA  M B  =  Ra  hs  Me  ,  (3-8) 

but  we  also  know  that  [v\a  =  g\  [v\c-  In  that  case  we  conclude 

A4 A  =  h'A  h b  '  (3-9) 

which  shows  how  orientation  matrices  chain  together.  The  process  can  be  continued  in¬ 
definitely  of  course:  g1^  =  g^  g^  g^,  and  so  on. 

Appendix  A  discusses  an  application  of  orientation  matrices  to  quantifying  the  various 
terms  used  in  discussions  of  wind  motion  over  a  ship. 

The  “Second  Fundamental  Question”  belongs  under  the  umbrella  of  how  to  rotate 
vectors,  which  we  look  at  next. 

3.2  Rotation  Mathematics 

Second  Fundamental  Question:  Suppose  S'  was  originally  coincident  with  S,  but 
now  has  been  orientated  according  to  g^ .  S'  took  a  snapshot  of  a  vector  u  and  carried  it 
along  to  make  a  new  vector  v.  We  ask:  how  is  [v]s  related  to  [ttjg?  First,  by  construction, 
the  components  of  v  in  S'  are  just  the  components  of  u  in  S: 


Ms'^Ms-  (3-10) 

In  that  case  we  can  immediately  write 

Ms  =  Us  Ms'  =  Us  Ms  >  (3-11) 

so  that  [v]s  relates  to  [m]^  via  the  orientation  matrix. 

3.2.1  How  to  Rotate  a  Vector  using  a  Matrix 

Now  refer  to  Figure  1,  in  which  a  vector  u  (or  rather,  a  snapshot  of  u )  is  rotated  in  a 
right-handed  sense  about  an  axis  vector  n  to  produce  v.  This  rotation  is  independent  of 
any  coordinate  system.  Write  it  as  Ren\ 

v  =  Renu.  (3.12) 
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Working  with  vectors  generally  requires  their  components  to  be  specified,  and  calculating 
the  components  of  a  rotated  vector  requires  a  coordinate  system.  The  core  result  of  rota¬ 
tion  analysis  is  that  in  S,  the  rotation  (3.12)  can  be  performed  by  a  matrix  multiplication 
using  the  matrix  representative  [R^]s  of  the  rotation: 

Ns  =  [Rn]s  Ns  >  (3.13) 

where  the  rotation  matrix  is,  for  unit-length  n, 

[Ren]s  =  1  +  sin0[n]*  +  (1  -cos0)  ([n]£)2.  (3.14) 


Here  the  first  “1”  is  the  3x3  identity  matrix.  The  cross  matrix  [n]g  appears  widely  in 
rotation  theory,  and  implements  the  vector  cross  product.  For  any  two  vectors  a,  b,  we 
have  [a]g  [b\^  =  [a  x  b]s,  where 
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(3.15) 


The  rotation  matrix  (3.14)  is  well  known  and  can  be  derived  by  combining  equations  (4.22) 
and  (4.25)  of  [4],  (There,  the  “[  ]g”  notation  of  this  report  is  suppressed  for  simplicity 
because  only  one  coordinate  system  is  used.)  See  also  [5],  whose  (1)  is  equivalent  to 
our  (3.14). 

A  rotation  can  be  inverted  (undone)  with  a  reverse  turn,  or  by  reversing  the  axis: 

(O"1  =R~d  =  Re_n,  (3.16) 


and  it  can  be  shown  easily  from  (3.14)  that  rotation  matrices  are  orthogonal,  meaning 


For  brevity,  we  write  a  rotation  about  an  axis,  say  R®  ,  as  Rex. 
A  rotation  R1  followed  by  R2  produces 


v  =  R2(R1u)  =  ( R2Ri)u , 


(3.17) 


(3.18) 


so  that  the  result  is  the  same  as  rotating  with  a  single  matrix  R2R\. 

Rotations  around  the  coordinate  axes  are  particularly  simple  and  known  as  Euler 
matrices: 
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El  =  [R*]s  = 
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(3.19) 


We  use  subscripts  1,  2,  3  for  the  Euler  matrices  because  they  depend  not  on  any  particular 
basis,  but  only  on  the  numerical  ordering  of  the  relevant  axes  within  the  basis  set.  That  is, 

Ef  =  [Rl]s  =  =  [R^js"  =  •  •  •  (3.20) 
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Given  two  sets  of  axes  S,S',  the  orientation  matrix  of  one  relative  to  the  other  is 
precisely  a  rotation  matrix: 

=[Ren\s-  (3-21) 

(We  use  this  identity  later  when  dealing  with  direction-cosine  representations.)  Prove  this 
by  starting  with  the  fact  that  /ij  gives  the  orientation  of  the  axes  S'  relative  to  the  axes  S. 
Now  consider  forming  the  S'  axes  by  rotating  a  snapshot  of  S  by  R so  that  the  x'  axis 
is  the  result  of  rotating  the  x  axis,  and  similarly  for  y  and  z.  The  orientation  matrix  is 


qi 


(3.22) 


Matrix  multiplication  proceeds  column  by  column — as  explained  ahead  in  (5.12) — and 
each  column  of  /i|  was  constructed  by  rotating  the  corresponding  coordinate  vector  of  S. 
That  is,  because 

Gx/  Rn  ^ x  ;  ^ y '  Rn  ^ y  ?  ^z'  ;  (3.23) 

we  can  immediately  write 


identity  matrix! 


(3.24) 


So  the  matrix  /i(f  that  gives  the  orientation  of  the  S'  axes  in  S  coordinates  is  precisely 
the  rotation  matrix  that  produced  the  S'  axes  from  the  S  axes.  This  is  an  important  fact 
that  aids  in  relating  orientation  and  rotation  matrices,  depending  on  which  is  more  easily 
constructed  in  the  problem  we  wish  to  solve. 

The  intricacies  and  richness  of  rotation/orientation  theory  arise  because  rotations 
around  different  axes  don’t  commute — meaning  they  cannot  be  swapped  without  changing 
the  result.  (They  do  commute  when  around  the  same  axis.)  The  extent  of  the  noncom- 
mutivity  is  given  by  examining  the  value  of 

[Rem,  Rt]  =RemRt~  Rt  RL  •  (3.25) 

If  that  value  is  zero,  the  rotations  commute,  but  in  general  [Rem,  Rn}  won’t  be  zero.  We 
can  calculate  its  value  using  (3.14),  omitting  the  “[  ]g”  notation  for  simplicity  since  this 
analysis  uses  a  single  coordinate  system.  Write 

[R9m,  Rn]  =  [l  +  sin#mx  +  (1  —  cos  8)mx2  ]  [l  +  sin cf>nx  +  (1  —  cos  0)nx'  ] 

—  [swap:  m  n,  6  <->  f] 

=  sin  6  sin</>  [mx,  nx]  +  (1  —  cos  6)  sin</>  [mx2,  nx] 

+  sin#  (1  —  cos  <f>)  [mx,  nx2  ]  +  (1  —  cos  9)(  1  —  cos^)  [mx2,  nx2  ] .  (3.26) 

As  6  and  f  become  small,  this  reduces  to 

\Rmi  Rn]  =  [m><> n>< ]  +  higher-order  terms  in  the  angles.  (3.27) 

That  is,  the  error  we  make  in  swapping  two  rotations  is  second  order  in  the  angles.  So  ro¬ 
tations  commute  to  first  order;  the  phrase  “infinitesimal  rotations  commute”  is  often  used, 
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meaning  that  we  can  swap  rotations  freely  when  working  with  infinitesimals.  (See  Sec¬ 
tion  6  for  more  on  this.)  Note  that  both  rotations  must  be  infinitesimal  if  we  wish  to  swap 
them. 

In  fact,  (3.27)  can  be  derived  more  quickly  than  we  have  just  done  in  the  last  paragraph. 
The  key  is  to  realise  that  (3.14)  can  be  written  in  the  compact  form 

K  =  e9nX,  (3.28) 

where  again  we  are  suppressing  the  “[  ]g”  notation.  This  exponential  is  defined  as  a  series 
in  the  usual  way.  If  we  now  calculate  [R6m,  Rn]  by  writing  each  of  its  rotations  as  an 
exponential  series,  then  the  final  small-angle  result  in  (3.27)  follows  in  a  straightforward 
way,  although  the  exact  expression  of  (3.26)  is  not  so  easily  obtained. 

According  to  (3.27),  whereas  rotations  through  small  angles  can  be  swapped  with  ever-increasing 
impunity  as  the  angles  are  made  still  smaller,  rotations  through  large  angles  cannot  be  swapped 
without  incurring  a  non- negligible  error.  We  might  attempt  to  swap  two  large  rotations  while  avoid¬ 
ing  the  above  error  of  order  6(f>,  by  dividing  each  into  infinitely  many  infinitesimal  rotations,  then 
swapping  adjacent  rotations  to  “bubble”  one  set  through  the  other.  In  fact,  this  procedure  fails. 

To  see  why,  start  with  two  large  rotations  through  9  and  (f>  around  distinct  axes,  and  divide 
each  into  N  tiny  rotations  8/N  and  <f>/N  (where  N  is  large).  Now  perform  the  required  N2  swaps, 
where  each  swap  adds  an  error  ~  9(j>/N2.  Thus  the  final  error  is  of  order  9(f),  and  we  have  gained 
nothing.  The  same  argument  shows  that  a  large  rotation  cannot  even  be  swapped  with  an  in¬ 
finitesimal  rotation.  The  PRT’s  pseudo-swap  of  large  rotations  works  only  because  it  changes  the 
rotation  axes. 


3.2.2  How  to  Rotate  a  Vector  using  a  Quaternion 

Euler  and  Rodrigues  are  credited  with  having  invented  a  set  of  four  numbers  that  can  be 
used  to  rotate  a  vector.  These  numbers  turn  out  to  be  identical  to  a  set  discovered  later 
by  Hamilton,  known  as  unit-length  quaternions.  We  will  briefly  cover  how  quaternions  are 
used;  the  relevant  proofs  can  be  found  in  [4].  A  quaternion  can  be  written  as  a  4-element 
vector,  typically  denoted  (a0,a),  where  the  last  three  elements  are  naturally  represented 
by  a  3-element  vector.3  In  fact,  quaternions  predate  vectors,  being  invented  in  the  mid 
19th  century  by  Hamilton  as  he  sought  to  generalise  the  idea  of  complex  numbers  a  +  bi. 
He  eventually  obtained  the  expression  a0  +  cqi.  +  a2j  +  o3k,  where  i,  j,k  had  the  neces¬ 
sary  algebraic  properties  to  allow  this  generalisation  to  use  a  minimum  of  new  ideas.4 
This  notation  eventually  split  into  a  “scalar”  part — the  first  component — and  a  “vector” 
part — the  last  three  components,  and  this  split  gave  rise  to  the  modern  notation  i,j,k 
for  generic  cartesian  basis  vectors.  So  Hamilton’s  a0  +  a\i  +  a2j  +  a3k  has  become  the 
modern  (a0,a).  We  must  use  more  descriptive  names  than  i,j,k  for  our  basis  vectors  in 
order  to  accommodate  more  than  one  coordinate  system:  eXl^Y^ez  f°r  ECEF,  ex,ey,  ez 
for  an  aircraft,  and  so  on. 

3We  are  using  the  symbol  for  a  proper  vector  a  to  represent  the  coordinate  vector  that  the  quaternion’s 
last  three  elements  can  be  defined  to  be.  The  distinction  is  understood  but  not  worth  dwelling  on  in  this 
context,  because  our  discussion  of  the  basic  properties  of  quaternions  uses  only  one  set  of  coordinates. 
After  all,  writing  our  quaternions  using  a  pre-defined  S  as  (a0,  [a]g) — or  really  (a0,  [a] 5) — would  just  be 
tedious. 

4For  Hamilton,  the  i  in  this  expression  generalised  the  complex  number  i,  but  from  the  point  of  view 
of  using  quaternions  for  rotation,  it’s  actually  the  quaternion  k  that  generalises  the  complex  number  i, 
because  both  of  these  implement  a  rotation  in  the  xy  plane. 
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Addition  and  scalar  multiplication  of  quaternions  is  defined  element  by  element,  just 
as  for  coordinate  vectors.  Quaternion  multiplication  is  associative  and  defined  by 

(o0,  a)  (b0,  b)  =  ( a0b0  —  a  b,  a0b  +  b0a  +  a  x  b) .  (3.29) 

The  quaternion  for  a  rotation  around  unit  vector  n  through  angle  6  is 

Qn  =  (cos  0/2,  n  sin  0/2) .  (3.30) 

It  can  be  shown  [4]  that  a  vector  u  is  rotated  through  9  around  unit  vector  n  to  produce 
vector  v  via  a  double  quaternion  multiplication: 

(0,v)  =  Qen(0,u)Q~e.  (3.31) 

[We  see  in  (3.31)  the  beginnings  of  just  why  it  was  that  the  last  three  elements  of  quater¬ 
nions  were  eventually  split  off  to  form  a  new  entity  called  a  vector:  because  the  quaternions 
that  are  used  in  place  of  vectors  u  and  v  are  (0,  u)  and  (0,  u).] 

Why  are  there  two  multiplications  in  (3.31)?  They  are  the  price  to  be  paid  for  the 
simple  structure  of  Qeni  versus  the  more  complicated  structure  of  the  matrix  Ren  in  (3.14). 
Notice  that  the  elements  of  n  occur  only  singly  in  Q^',  contrast  this  with  Ren ,  where  the 
last  term  in  (3.14)  has  products  of  the  elements  of  n  sitting  inside  nxL  Quaternions’ 
sparseness  gives  them  simplicity,  but  the  price  to  be  paid  is  that  two  multiplications  are 
needed  to  incorporate  those  required  products  of  the  elements  of  n,  in  order  to  achieve  the 
same  as  a  matrix  multiplication  by  R^-  This  double  multiplication  is  also  why  (3.30)  uses 
just  half  the  rotation  angle. 

A  quaternion  is  defined  to  have  length 

|(a0,a)|  =  y'ao  +  |a|2' ,  (3.32) 

in  which  case  it  follows  that  all  rotation  quaternions — that  is,  of  the  sort  (3.30) — have  unit 
length.  This  corresponds  to  the  unit  determinant  of  rotation  matrices. 

The  identity  quaternion  for  rotation — corresponding  to  no  rotation  at  all — is  Q®n  =  (1,0), 
being  a  rotation  through  zero  angle  around  any  axis.  Because  (3.30)  makes  it  trivial  to 
build  the  quaternion  associated  with  an  angle  and  an  axis,  quaternions  have  traditionally 
come  to  be  treated  as  the  way  to  implement  the  angle-axis  approach  to  quantifying  an 
orientation.  But  we  should  not  forget  that  all  they  do  is  rotate  a  vector,  just  as  does  the 
matrix  in  (3.14). 

The  inverse  of  Q ^  is 

{Q^y1  =  Qe-n  =  Q-e,  (3.33) 

and,  just  as  for  matrices,  the  product  of  quaternions  Qi  and  Q2  satisfies 

{Q1Q2)  =  Q2  Qi  ■  (3.34) 

Equations  (3.31)-(3.34)  show  how  successive  rotations  are  written  with  quaternions:  a  ro¬ 
tation  Qi  followed  by  Q2  produces 

(0,  v)  =  Q2Q\  (0,  u)  Qx  lQ2  1  =  Q2Qi  (0,  u )  (Q2Q\)  1 ,  (3.35) 

so,  just  as  in  (3.18),  the  result  is  the  same  as  rotating  with  a  single  quaternion  Q2Q\. 
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3.2.3  The  Rotation  Matrix  and  Rotation  Quaternion  in  Quantum 
Mechanics 

The  exponential  form  (3.28)  of  the  rotation  matrix  is  derived  in  some  quantum  mechanics 
textbooks,  because  it’s  used  in  that  subject  to  deal  with  angular  momentum,  which  is  a 
property  of  fundamental  particles.  If  you  do  refer  to  a  quantum  mechanics  text,  you’ll  find 
a  different  notation,  perhaps  accompanied  by  an  obfuscated  way  of  deriving  (3.28).  We’ll 
pay  a  brief  visit  to  that  notation  here  to  demystify  it. 

First,  omit  the  explicit  [  ]g  on  all  vectors  for  simplicity,  and  use  the  linearity  of  the 
cross  operator  to  expand  the  cross  matrix  into  three  component  matrices:5 

Tl  (jlx€ix  T  TlySy  T  TlZGZ  )  TlX£X  T  TlySy  T  *  (3.37) 

In  quantum  mechanics  the  matrices  ex  ,  e*  ,  ez  are  denoted  —iJx,  — iJy,  — i Jz  respectively, 
and  the  matrix  triplet  (Jx,  Jyl  Jz)  is  written  as  J,  as  though  it  were  a  vector.  Hence  nx  is 
written  —in  J,  and  the  rotation  matrix  becomes 

Ren  =  eenX=e~WnJ.  (3.38) 

The  only  reason  for  the  introduction  of  I  is  that  the  matrices  Jx,Jy,  Jz  eventually  become 
identified  with  angular  momentum  in  a  quantum  mechanical  sense. 

It’s  notationally  efficient  to  derive  (3.28)  as  e0nX ,  and  only  then  to  branch  off  into 
quantum  mechanics  through  the  use  of  the  angular  momentum  matrices.  Instead  of  doing 
that,  the  standard  quantum  mechanical  treatment  introduces  Jx,Jy,  Jz  at  the  very  start  of 
the  subject-  -  before  the  rotation  matrix  has  even  been  derived—  by  writing  ex  as  —  iJx  and 
so  on;  it  thus  presents  a  real,  geometrical,  analysis  using  complex  numbers  that  all  cancel 
internally  (since  the  J  matrices  are  pure  imaginary).  The  effect  is  rather  like  writing 
2  x  3  =  21  x  — 3i  =  6.  There  is  no  pedagogical  value  in  such  a  procedure.  Of  course, 
quantum  mechanics  requires  that  Jx,Jy ,  Jz  involve  complex  numbers;  that  is  not  the  issue.6 
Rather,  why  introduce  I  into  a  rotational  analysis  where  it’s  not  needed?  There  is  nothing 
deep  or  elegant  about  writing  a  rotation  as  e~l9n  J  before  introducing  the  concept  of 
angular  momentum,  and  I  suspect  that  most  users  of  the  quantum  mechanical  notation 
would  be  hard  pressed  to  actually  rotate  a  vector  using  e~xdn' J — assuming  they  were  aware 
that  this  is  what  e~x6n  J  actually  does — despite  the  pages  of  obscure  analysis  that  you’ll 
find  on  the  subject  in  quantum  mechanics  textbooks.  If  you  do  refer  to  the  quantum 
mechanics  approach  for  some  rotation  theory,  recognise  this  unnecessary  comp  lex- number 
baggage  for  what  it  is,  and  the  analyses  will  become  a  little  less  confusing. 

Similar  comments  apply  to  the  quaternions  used  in  quantum  mechanics.  It  turns  out 
that  a  quaternion  can  be  represented  by  a  2  x  2  complex  matrix  [4] ,  and  you  will  find  them 

5Equation  (3.37)  shows  the  compactness  of  the  cross  notation.  We  could  also  have  used  its  full  matrix 
form  to  get  the  same  result: 

nz  ny 
0  —nx 
nx  0 


6  On  this  note,  given  that  quantum  mechanics  is  based  on  the  complex  Fourier  transform,  and  given 
that  the  complex  Fourier  transform  can  always  be  replaced  by  the  real  Hartley  transform  [6],  perhaps  the 
last  word  has  yet  to  be  said  on  this  subject.  But  that  doesn’t  affect  my  comments  above. 
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written  that  way  in  quantum  mechanics  where  the  fundamental  matrices  used,  analogous 
to  the  basic  quaternions  l,i, j,k,  are  called  Pauli  matrices.  Just  as  for  the  Jx,Jy,Jz 
matrices  in  the  previous  paragraphs,  Pauli  matrices  are  highly  useful  in  describing  angular 
momentum.  But  being  complex,  they  need  to  be  associated  with  an  extra  factor  of  i  (the 
square  root  of  —1,  not  the  quaternion)  in  order  to  be  used  to  rotate  a  vector,  so  it’s  best 
to  avoid  this  quantum  mechanics  notation  if  you  want  to  learn  rotation  theory. 

3.2.4  Rotating  One  Point  About  Another 

Finally,  consider  rotating  the  point  (10,0,0)  by  90°  about  an  axis  (0,0,1)  that  passes 
through  the  point  (9,0,0).  The  result  is  the  point  (9, 1,0),  which  is  easily  seen  from  the 
simple  geometry  involved.  How  might  we  use  the  theory  of  rotation  to  rotate  points  in 
general? 

To  see  how,  we  should  be  precise  in  describing  the  foregoing  rotation  using  vector 
language.  We  rotated  the  point  described  by  the  position  vector  (10,0,0)  (necessarily 
quantified  as  a  coordinate  vector ,  a  notion  that  we’ll  take  as  given!)  by  first  forming  a 
proper  vector :  the  arrow  from  (9,0,0)  to  (10,0,0).  This  arrow,  (1,0,0),  was  then  rotated 
by  90°  about  the  axis  proper  vector  (0,0,1)  to  give  (0,1,0).  Finally  this  arrow  (proper 
vector)  was  added  to  the  position  vector  (9, 0,  0)  to  give  the  position  vector  for  the  sought- 
after  point  (9, 1, 0). 

The  general  procedure  is  straightforward  to  write  down.  We  require  to  rotate  the 
point  x  by  angle  8  about  the  point  p,  through  which  runs  the  axis  n  (a  unit  vector).  The 
resulting  sought-after  point  is  denoted  x' .  We  simply  rotate  the  proper  vector  x  —  p  to 
produce  the  proper  vector  x'  —  p: 

x'  -  p  =  Ren  (x  -  p) ,  (3.39) 

so  that  the  required  point  is  x'  =  p  +  Ren  (x  —  p). 

4  Quantifying  Orientation 

We’ve  seen  that  an  object’s  orientation  can  only  be  expressed  relative  to  some  given  base 
orientation,  and  is  conveniently  written  within  some  chosen  coordinate  system  S.  It  can 
be  quantified  by  embedding  a  set  of  axes  within  the  object,  calling  these  the  S'  axes,  and 
specifying  the  coordinates  of  the  S'  basis  vectors  in  S.  These  coordinate  vectors  are  written 
as  the  columns  of  the  orientation  matrix  ,  also  called  the  direction-cosine  matrix. 

But  in  (3.24)  we  saw  that  p,g'  is  in  fact  a  rotation  matrix:  we  can  interpret  that  equation 
as  saying  that  rotates  snapshots  of  the  S  axes  to  produce  the  S'  axes.  This  allows  for 
alternative  methods  of  quantifying  an  object’s  orientation,  by  constructing  the  orientation 
through  imagining  that  the  object’s  axes  were  originally  coincident  with  those  of  S,  and 
we  rotated  it  in  some  chosen  way  to  arrive  at  its  final  orientation.  The  three  methods  that 
are  used  commonly  are  described  below.  It’s  important  to  understand  that  they — and  any 
others  that  could  be  added — all  ultimately  do  the  same  thing :  they  construct  the  S'  basis 
vectors  from  the  S  basis  vectors.  The  choice  of  one  method  over  another  is  partly  governed 
by  the  parameters  of  the  problem  being  solved,  and  partly  by  the  nuances  of  any  necessary 
computer  coding. 
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Orientation  or  Direction-cosine  matrix:  As  noted  after  (3.4),  these  are  alternative 
names  for  the  matrix  whose  columns  are  the  basis  vectors  of  S'  coordinated  in  S. 
This  matrix  completely  specifies  the  object’s  orientation. 

Angle-axis:  Here  we  specify  an  angle  6  and  a  unit-length  axis  vector  n.  rotating  snap¬ 
shots  of  the  S  basis  vectors  through  9  around  n  produces  the  S'  basis  vectors. 

Euler  angles:  We  specify  the  angles  of  a  sequence  of  rotations  that  rotates  a  snapshot 
of  the  S  axes  onto  the  S'  axes;  typically  3  rotations  are  used  around  the  coordinate 
axes;  these  angles  are  then  called  Euler  angles.  The  need  for  three  rotations  rather 
than  the  angle-axis  method’s  single  rotation  is  due  to  the  fact  that  the  Euler  angles 
specify  rotations  around  a  prescribed  set  of  cartesian  axes.  In  contrast,  the  angle- 
axis  method  manages  with  just  one  axis — but  that  axis  will  seldom  coincide  with 
any  “simple”  axis  such  as  the  cartesian  axes. 

The  first  method  is  specifically  a  matrix  method.  The  orientation  matrix  is  also  a  rotation 
matrix,  as  we  saw  in  (3.24).  The  last  two  methods  involve  rotation,  and  a  rotation  can  be 
done  using  either  a  matrix  or  a  quaternion.  The  choice  of  which  of  these  two  tools  to  use  is 
ours;  the  methods  stand  above  the  tools  used  to  implement  them.  But  the  important  point 
to  remember  is  that  the  methods  all  do  the  same  thing,  and  in  particular  the  matrices  that 
can  be  produced  from  all  three  of  them  will  be  identical  (up  to  numerical  inaccuracies). 

Comparing  methods:  Given  that  the  orientation  matrix  is  identical  to  using  a  single 
rotation  to  describe  an  orientation,  the  choice  of  method  generally  comes  down  to 
using  either  one  or  three  rotations  to  rotate  the  base  orientation.  Using  a  single 
rotation  is  as  lean  and  simple  as  things  can  be.  It  presents  no  great  numerical 
difficulties  and  has  now  become  a  proven  approach.  Using  three  Euler  angles  does 
have  its  difficulties:  the  rotations’  lack  of  commutivity  means  that  users  of  Euler 
angles  often  find  them  confusing,  and  there  are  numerical  problems  with  using  three 
angles  that  we’ll  investigate  in  Section  5.4.  Nevertheless,  rotating  by  Euler  angles 
corresponds  to  the  yaw-pitch-roll  language  and  manoeuvres  familiar  to  pilots,  and 
this  means  that  there  will  always  be  a  place  for  Euler  angles  in  orientation  theory, 
at  least  when  applied  to  flight. 

Comparing  tools:  The  price  we  pay  for  the  simplicity  of  matrix  multiplication  is  that 
rotation  matrices  have  redundancy  in  their  elements.  A  quaternion  can  be  viewed  as 
a  way  of  “rendering  down”  a  rotation  matrix  to  its  core  information,  but  this  comes 
with  a  more  complicated  form  of  multiplication  than  that  of  matrices. 

There  is  an  important  computational  point  to  be  considered  here:  if  we  must 
extrapolate  orientation  information  (as  in  Section  6),  we’ll  be  updating  matrices 
or  quaternions  using  numerical  algorithms.  Throughout  such  procedures,  we  must 
always  ensure  that  the  matrices  produced  have  determinant  1,  as  all  rotation  matrices 
do,  or  that  the  quaternions  produced  have  length  1,  as  all  rotation  quaternions  do. 
This  mitigates  problems  due  to  numerical  round-off  error;  if  we  didn’t  at  least  ensure 
our  matrices  or  quaternions  had  this  property,  then  the  results  of  rotating  vectors 
would  soon  begin  to  look  skewed. 

A  suitable  replacement  for  a  quaternion  Q  whose  length  has  strayed  from  1 
is  to  replace  Q  by  itself  divided  by  (meaning  each  of  its  elements  by)  its  length, 
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calculated  in  (3.32):  Q~^>Q/\Q\-  This  is  a  simple  numerical  exercise.  In  con¬ 
trast,  re-orthogonalising  a  rotation  matrix  R  requires  significantly  more  compu¬ 
tation.  A  suitable  definition  for  how  “close”  its  replacement  R'  should  be  leads 
to  R'  =  R{RtR)^1/2 ,  which  can  be  evaluated  using  singular  value  decomposition. 
That  is,  first  use  that  technique  to  write  R  =  U SVl  for  matrices  U,  S,  V  as  outlined 
in  singular  value  decomposition  theory;  their  properties  quickly  produce  the  simple 
expression  R'  =  UV1. 

There  persists  a  view  in  the  orientation  community  that  angle-axis  goes  hand  in  hand  with 
quaternions,  whereas  Euler  angles  are  involved  with  matrices.  Because  of  this,  the  tool 
(quaternions)  used  to  implement  the  angle-axis  description  of  an  orientation  has  gradually 
come  to  be  compared  with  the  Euler  angles  method  (three  rotations).  But  that’s  comparing 
apples  and  oranges;  the  three  methods  described  above  should  be  compared  on  one  level, 
and  the  two  tools  (matrices  and  quaternions)  should  be  compared  on  another  level. 


5  Alternative  Approaches  to  the  Applications  of 

Reference  [1] 

Sections  4. 1-4.3  of  [1]  gave  examples  of  using  rotations  to  construct  orientation  scenarios. 
We’ll  now  rework  those  calculations  using  the  notation  of  this  report.  The  result  will  be 
a  more  streamlined  approach  to  the  problems  being  solved. 

5.1  Rotation  and  Flight  over  a  Sphere 

Any  translation  that  a  body’s  centroid  undergoes  while  the  body’s  orientation  is  changing 
does  not  affect  that  orientation.  An  aircraft  that  rolls  to  fly  inverted  is  inverted,  irre¬ 
spective  of  the  fact  that  it’s  moving  sideways  as  it  turns  upside  down.  We  can  use  this 
fact  to  our  advantage  when  visualising  rotation.  For  example,  picture  how  an  aircraft’s 
orientation  changes  as  it  pitches  nose  down,  rotating  through  —90°  about  its  starboard 
wing.  We  can  arrive  at  the  same  final  orientation  by  imagining  the  aircraft  to  be  anywhere 
on  a  sphere — say,  a  point  on  the  sphere’s  equator.  Now  allow  the  aircraft  to  move  over 
the  sphere  (like  a  car  being  driven,  always  touching  the  sphere  so  that  it’s  straight  and 
level  at  every  point)  until  it  has  flown  to  the  North  Pole.  This  change  in  orientation  is 
precisely  a  rotation  through  —90°  around  its  starboard  wing;  the  motion  over  the  sphere 
is  of  no  consequence.  This  is  a  simple  example,  but  we  can  picture  a  sequence  of  more 
complicated  rotations  in  the  same  way,  by  tracing  the  aircraft’s  path  over  the  sphere. 

5.2  Constructing  Local  North-East— Down  Coordinates 

Section  4.1  of  [1]  posed  the  question  “If  we  are  in  Adelaide  and  Earth  is  transparent,  what 
is  the  compass  bearing  of  Brussels  if  we  are  looking  straight  at  it  through  Earth?”.  We 
can  answer  this  question  by  finding  the  local  north  and  east  components  of  the  vector  rBy^ 
pointing  from  Adelaide  A  to  Brussels  B,  and  using  these  with  some  simple  trigonometry 
to  find  the  required  angle.  Calculate  these  components  by  constructing  the  local  north¬ 
east-down  (NED)  axes  at  Adelaide;  this  is  done  by  calculating  the  coordinates  of  those 
axes’  basis  vectors  in  Earth-centred,  Earth-fixed  coordinates  (ECEF),  defined  shortly. 
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Figure  2:  The  local  north-east-down  vectors  at  a  given  point  on  Earth’s 
surface  can  be  constructed  by  rotating  the  ordered  set  of  vectors  ez,eY,  —ex 
(at  zero  lat/long,  lower  left)  to  produce  the  ordered  set  n,e,d  (upper  right) 


First,  we’ll  construct  the  vector  rBj 4  pointing  from  Adelaide  A  to  Brussels  B,  specifying 


its  components  in  ECEF  coordinates  E.  ECEF  coordinates  are  cartesian  with  origin  at 
Earth’s  centre,  using  the  X,Y,  Z  axes  of  Fig.  2.  They  often  simplify  geodesy  calculations. 
The  ECEF  system  has  the  following  unit-length  basis  vectors: 

1.  ex  points  from  Earth’s  centre  through  the  Equator  at  a  =  uj  =  0°, 

2.  eY  points  from  Earth’s  centre  through  the  Equator  at  a  =  0°,lu  =  90°,  and 

3.  ez  points  from  Earth’s  centre  through  the  North  Pole. 

The  proper  vector  rBj 4  pointing  from  Adelaide  to  Brussels  is  found  from  Brussels’ 

ECEF  position  vector  rB(Z  (which  is  relative  to  Earth’s  centre  C )  and  Adelaide’s  ECEF 
position  vector  rAC: 


(5.1) 


rBA  ~  rBC  ~  rAC 


where  each  of  these  vectors  is  calculated  in  the  following  way.  Given  a  point  P’s  latitude  a, 


longitude  u,  and  height  h  above  the  WGS-84  ellipsoid,  its  ECEF  position,  [ rPC\E ,  is  the 
(. X ,  Y.  Z )  triplet  of  the  point  using  the  axes  convention  of  Fig.  2.  Using  Earth’s  semi-major 
and  semi-minor  axis  lengths  in  the  WGS-84  scheme, 


semi-major:  a  =  6,378,137  m, 
semi-minor:  b  =  6,356,752.3142  m, 


(5.2) 


an  analysis  of  the  shape  of  an  ellipsoid  yields  P’s  ECEF  coordinates  as 
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lb2  \ 

Z  =  I  |  +  h  1  sin  a  .  (5-3) 

\V«2  cos2  a  +  b2  sin2  a  ) 

(This  equation  is  equivalent  to  equation  (2.2)  of  [1],  but  is  written  here  in  a  way  that 
treats  a  and  b  more  symmetrically.)  Assume  the  two  cities  are  at  h  =  0  and  use 

Adelaide:  latitude  a  =  —34.9°,  longitude  u>  =  138.5°, 

Brussels:  latitude  a  =  50.8°,  longitude  u  =  4.3°. 

Equation  (5.3)  converts  these  parameters  into  two  position  vectors  of  the  cities  with  respect 
to  Earth’s  centre  C.  These  positions  of  Adelaide  and  Brussels  are,  respectively, 


'-3.922' 

'4.028' 

rAc\E  ~ 

3.470 

-3.629 

x  106  metres, 

[vbc\e  — 

0.303 

4.920 

x  106  metres, 


(5.4) 


and  a  subtraction  gives  the  vector  pointing  from  Adelaide  to  Brussels  (with  all  internal 
calculations  in  this  report  taken  to  a  large  number  of  decimal  places  and  then  rounded  to 
3  decimal  places  when  written  down): 


\Tba\e  ~  \ rBC  ~  r Ac\e  ~  \rBc\E  ~  [r Ac\e  ~ 


7.950 

-3.167 

8.548 


x  10b  metres. 


(5.5) 


We  have  [rBj^\E,  but  require  the  first  two  components  of  [rB^]N,  since  these  are  the  north 
and  east  components  needed  to  calculate  the  required  bearing  of  Brussels  as  seen  from 
Adelaide.  These  two  coordinate  vectors  are  related  by 


[vba]n  ~ — -  [tba\e  ■  (5-6) 

The  ECEF  is  particularly  easy  to  calculate  within,  so  we  might  calculate  /r®  by  first 
finding  //)) .  For  this,  we  need  the  NED  coordinate  system  N  at  a  given  point  on  Earth  at 
latitude  a  and  longitude  <n,  also  shown  in  Figure  2.  The  NED  coordinate  axes  N  are  defined 
by  unit  basis  vectors  n  (pointing  north),  e  (pointing  east),  d  (pointing  down,  =  n  x  e). 
We  require  to  find  these  basis  vectors,  expressed  in  ECEF  coordinates.  That  is,  we  require 
Wei  [e].E;  [d\E- 

The  flying-aircraft  analogy  of  Section  5.1  can  be  used  here.  Earth  is  approximately  an 
oblate  spheroid — not  exactly  a  sphere;  but  the  geodetic  latitude  a  in  Figure  2  is  defined 
in  such  a  way  that  we  can  treat  Earth  as  exactly  spherical  when  using  a  to  construct 
NED  axes.  The  reason  is  because  the  “down”  vector  at  any  latitude  a  is  defined  to  be  the 
result  of  rotating  the  “down”  vector  at  the  Equator  through  the  angle  a  (around  local  west , 
which  allows  latitude  to  increase  towards  the  North  Pole,  purely  by  historical  convention). 
This  is  the  same  resulting  orientation  of  the  “down”  vector  that  would  result  by  rotating  it 
simply  by  an  angle  a — or  by  flying  an  aircraft  through  an  arc  of  angle  a  in  a  great  circle 
around  a  spherical  Earth. 

The  vectors  n,e,d  are  constructed  by  rotating  eZleY,—ex  respectively  twice,  corre¬ 
sponding  to  implementing  the  changes  in  latitude  and  longitude: 


{eZi  eYi  ~ex}  -►  in,e,d}  . 


(5.7) 
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Figure  3:  Two  methods  for  implementing  the  vector  construction  of  (5.7). 
We  can  rotate  the  original  ordered  set  {ez,eY,—ex}  first  around  —eY 
through  a  (or,  equivalently,  around  eY  through  —a)  to  visit  point  P,  and 
then  around  ez  through  u.  Alternatively,  we  can  rotate  the  original  ordered 
set  {ez,eY,—ex}  first  around  ez  through  ui  to  visit  Q,  and  then  around 
the  rotated  snapshot  of  eY  through  —a,  and  then  use  the  Pseudo- Reversing 
Theorem  to  avoid  having  to  calculate  the  rotated  snapshot  of  eY . 


Working  in  coordinates  E,  the  relevant  coordinate  vectors  are 

|  |[n]E,  [e]s,  .  (5.8) 

Two  ways  of  doing  the  rotations  are  immediately  apparent  in  Figure  2.  Either  of  these 
can  be  used,  and  we’ll  do  the  calculation  for  each  way  to  contrast  them. 

First  Way  to  Implement  (5.7):  Recalling  the  positive  conventions  for  latitude  and 
longitude,  rotate  each  vector  of  the  set  {ez,  eY,  —ex}  first  around  —  eY  through  a  (or, 
equivalently,  around  eY  through  —a)  to  arrive  at  point  P  in  Figure  3,  and  then  around  ez 
through  uj  to  become  its  corresponding  vector  in  {n,e,d}.  The  relevant  rotations  are 

ieZi  eYi  ~ ex}  Rya  *  R-z  (n> e>  ■  (^-9) 

Second  Way  to  Implement  (5.7):  Alternatively,  we  can  rotate  each  initial  vector  first 
around  ez  through  u  to  arrive  at  point  Q  in  Figure  3,  and  then  around  the  rotated  snapshot 
of  eY  through  —a: 

{eZi  eYi  ~  exj  ~ Rz  R(y)  in’  e’ ■  (5.10) 

But  the  Pseudo-Reversing  Theorem  allows  us  to  avoid  having  to  calculate  the  rotated 
snapshot  of  eY\  That  is,  applying  it  to  the  rotations  in  (5.10)  produces  (5.9)  straightaway. 
You  might  say  that  there  was  no  need  to  use  the  Pseudo-Reversing  Theorem  here,  since  we 
could  always  have  chosen  the  First  Way  above  of  doing  the  rotations.  That’s  true  for  this 
simple  example,  but  when  more  than  just  two  rotations  are  being  used,  the  Second  Way  is 


r 

0 

0 

-1 

0 

1 

0 

l 
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0 
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usually  the  most  natural  way  of  approaching  the  rotations — and  almost  certainly  the  only 
easily  visualised  way;  and  then  the  Pseudo-Reversing  Theorem  becomes  indispensable  to 
avoid  having  to  calculate  the  many  intermediate  snapshots  that  would  be  required  if  we 
were  to  implement  the  Second  Way  directly. 

The  Second  Way  actually  implements  a  procedure  from  differential  geometry  called  parallel  trans¬ 
port :  when  we  slide  the  north  and  east  vectors  from  zero  lat/long  via  point  Q  in  Figure  3,  we 
are  parallel  transporting  them;  the  notion  doesn’t  apply  to  the  down  vector,  but  it’s  always  the 
cross  product  of  north  and  east  anyway.  Parallel  transport  is  all  about  movement  along  geodesics, 
which  are  the  great  circles  on  a  sphere.  Because  great  circles  are  the  tracks  we  would  fly  (or  drive) 
if  we  didn’t  turn  a  steering  wheel,  they  are  natural  tracks  on  a  surface,  and  this  almost  certainly 
accounts  for  their  ease  of  use  in  visualising  rotations. 


The  two  rotations  of  (5.9)  now  allow  us  to  implement  (5.8)  using  matrices.  The  rele¬ 
vant  matrices  are  just  rotations  about  the  axes,  which  are,  of  course,  the  Euler  matrices 
of  (3.19): 

|  E-°  -  E%  -  {[n]E,  [e]E,  [ d]E }  .  (5.11) 

The  three  vectors  can  be  multiplied  together,  somewhat  in  parallel.  This  is  because  matrix 
multiplication  proceeds  column  by  column.  For  suppose  we  have  a  matrix  A  for  which 


f 

0 

0 

-1 

0 
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0 
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0 

n 

'1' 
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'10' 
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A 

2 

= 

b 
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20 

= 

q 

3 

c 

30 

r 

Then  we  can  immediately  write 


'1  10' 

a  p 

A 

2  20 

= 

b  q 

1 

CO 

CO 

o 

1 _ 

c  r 

(5.12) 


(5.13) 


This  “block  multiplication”  allows  us  to  condense  the  three  multiplications  of  coordinate 
vectors  into  one  matrix  multiplication  by  coalescing  each  set  of  three  coordinate  vectors 
into  a  matrix: 


[n\E  [e}E  [d]E 


_  rpoi  jp-at 

~  ^3  ^2 

COS  LO 
=  sin  lo 
0 


0  0-1 
0  1  0 
10  0 

—  sin  lo  0 
cos  lo  0 

0  1 


cos  a 

0 

—  sin  a 

'0 

0 

-1' 

0 

1 

0 

0 

1 

0 

sin  a 

0 

cos  a 

1 

0 

0 

—  cos  lo  sin  a 

—  sin  lo  sin  a 

cos  a 


—  sin  lo  —  cos  lo  cos  a 
cos  lo  —  sin  lo  cos  a 
0  —  sin  a 


(5.14) 


The  ECEF  coordinate  vectors  for  the  local  north,  east,  and  down  vectors  at  a  specific  loca¬ 
tion  on  Earth,  such  as  Adelaide,  are  now  easily  found.  This  city  has  latitude  a  =  —34.9°, 
longitude  lo  =  138.5°.  Use  these  numbers  to  write  (5.14)  as 


E e 


n 


Mi 


-0.429  -0.663  0.614 

0.379  -0.749  -0.543 
0.820  0  0.572 


(5.15) 
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so  that  the  ECEF  coordinates  of  the  unit  vector  pointing  to  local  north  at  Adelaide  are 
(—0.429,  0.379,  0.820),  and  similarly  for  east  and  down. 

Remember  now  that  the  transpose  of  ^  is  so  that  we  have  now  also  obtained  the 
ECEF’s  basis  vectors  in  NED  coordinates: 


iex]N  [ey]7V  iez]N  —  Mat  —  (Re)  —  [n]-E  [e]-E  M 


1 1 


(5.16) 


So,  for  example,  the  NED  coordinates  of  ex  are  (—0.429,  —0.663,  0.614).  Equation  (5.15) 
has  the  same  information  as  equation  (4.8)  in  [1], 

Alternative  methods  for  producing  the  orientation  matrix  (5.14)  from  rotations  exist, 
depending  on  what  vectors  we  take  as  our  initial  ordered  set  that  map  onto  n,e,d.  For 
example,  choosing  initial  vectors  ex,  eY,ez  to  map  respectively  to  n,  e,  d  alters  the  matrix 
multiplication  of  (5.14)  to  E%  ,  which  is  completely  consistent  with  the  rightmost 

matrix  in  the  first  line  of  (5.14)  being,  in  fact,  E^90  ■ 

Finally,  the  required  NED  coordinates  of  the  vector  pointing  from  Adelaide  to  Brus¬ 
sels  are 


\tba] 


N 


(3.7)  e  r_  n  (3.5) 


fJ-N  \VBa]e  = '=  {^eY  \rBA]E 


(5.15) 


-0.429  0.379  0.820 

-0.663  -0.749  0 

0.614  -0.543  0.572 


7.950' 

-3.167 

8.548 

x  106  metres 


2.403 

-2.896 

11.495 


x  106  metres. 


(5.17) 


The  first  two  components  of  this  vector  were  given  in  equation  (4.9)  of  [1],  and  we  can 
see  that  the  calculations  of  the  two  reports  agree.  The  bearing  of  Brussels  as  seen  in  a 
straight  line  through  a  transparent  Earth  is  given  by 


—  tan  1 


2.896 

2.403 


~  -50.3°  . 


(5.18) 


5.3  The  Direction  in  which  a  Pilot  Sights  an  Aircraft 

The  example  in  Section  4.2  of  [1]  placed  an  aircraft  over  Adelaide,  flying  north-east,  pitched 
up  at  20°  at  an  altitude  of  30,000  metres.  The  following  question  was  posed.  We  who 
pilot  the  aircraft  sight  another  aircraft  flying  over  Sydney  at  the  same  altitude.  In  which 
direction  in  our  cockpit  must  we  look  to  see  that  aircraft? 

This  calculation  is  very  similar  to  that  of  Section  5.2.  We  construct  the  vector  pointing 
from  us  to  the  other  aircraft,  then  produce  its  coordinates  using  our  aircraft’s  axes  as 
opposed  to  the  NED  axes  of  Section  5.2.  The  coordinate  vector  pointing  from  us  to  the 
other  aircraft  is  constructed  in  the  same  way  as  was  done  for  Adelaide  and  Brussels  in  the 
previous  section:  begin  with  the  locations  of  the  Adelaide  and  Sydney  aircraft: 

Adelaide:  latitude  a  =  —34.9°,  longitude  uj  =  138.5°,  height  30,000  m 
Sydney:  latitude  a  =  —33.9°,  longitude  ui  =  151.2°,  height  30,000  m. 
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y 


Figure  Conventional  labels  of  the  cartesian  body  axes  of  an  aircraft,  with 
directions  of  positive  rotation  indicated 


Now  apply  (5.3)  to  find  the  position  vector  of  the  Adelaide  aircraft  A  with  respect  to 
Earth’s  centre  C,  in  ECEF  coordinates  E,  and  likewise  for  the  Sydney  aircraft  S: 


'-3.941' 

'-4.666' 

vac\e  - 

3.486 

-3.646 

x  106  metres, 

irSc\E  ~ 

2.565 

-3.554 

x  106  metres.  (5.19) 


The  vector  pointing  from  us  in  the  Adelaide  aircraft  to  the  Sydney  aircraft,  in  ECEF 
coordinates  E,  is 


[rSA\E  -  [ vSc\e  ~  [rAC\E  — 


-7.252 

-9.213 

0.920 


x  10°  metres. 


(5.20) 


We  need  only  convert  this  vector  to  our  aircraft’s  coordinates.  Call  these  coordinates  uUs ”, 
so  that  we  require 

lrSA\us  =  hnjs  \tSa\e  •  (5.21) 

Given  (5.20),  it  remains  to  find 


t4s  = 


[' ex\  Us  [er]  Us  [ ez ]  Us 


(5.22) 


As  before,  label  the  ECEF’s  axes  X,  Y,  Z,  and  now  label  our  aircraft’s  axes  x,  y ,  z.  These 
are  shown  in  Figure  4.  There  are  two  main  ways  that  we  can  list  the  procedure  to  con¬ 
struct  n  :  one  is  easy  to  visualise  but  needs  the  Pseudo- Reversing  Theorem,  while  the 
other  is  perhaps  harder  to  visualise,  but  doesn’t  need  the  Pseudo-Reversing  Theorem. 
There  are  other  ways  that  permute  the  steps  we  outline,  but  they  are  just  more  of  the 
same  and  needn’t  be  considered  here. 


First  Way  to  Construct  /x^:  Work  in  the  ECEF,  so  invoke  g^s  =  to  allow  us 

first  to  calculate 


Us 


h'E 


(5.23) 
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This  might  be  more  naturally  constructed  than  y^js,  because  /rj?  refers  directly  to  the 
aircraft’s  axis  vectors  ex,ey,ez,  and  these  are  easily  constructed  by  rotating  snapshots 
of  the  ECEF’s  X,  Y,  Z  axes  to  “build”  our  aircraft’s  x,  y,  z  axes  in  the  following  way. 
Imagine  the  aircraft  initially  at  zero  lat/long,  with  its  x,y,z  axes  parallel  to  the  ECEF’s 
X,  Y,  Z  axes  respectively,  balancing  with  its  tail  on  the  ground  and  with  starboard  wing 
pointing  east.  In  the  following  rotations,  imagine  rotating  the  actual  aircraft  to  bring 
it  to  the  required  orientation  at  the  required  place  on  Earth — recalling  the  comments  of 
Section  5.1  and  the  discussion  of  geodetic  latitude  in  Section  5.2.  Start  with  snapshots 
of  the  ECEF’s  basis  vectors  ex,eY,ez  which  will  be  embedded  in  the  aircraft’s  body  to 
eventually  become  the  required  ex,ey,ez  vectors.  Now: 

1.  rotate  this  “base”  aircraft  through  —90°  around  its  starboard  wing,  being  the  ECEF’s 
Y  axis.  This  produces  a  new  set  of  “latest”  basis  vectors,  orientated  so  that  the  latest 
(rotated)  snapshot  of  ex  (embedded  in  the  aircraft’s  nose)  points  north,  the  latest 
(rotated)  snapshot  of  eY  (the  starboard  wing)  points  east,  and  the  latest  (rotated) 
snapshot  of  ez  points  down. 

2.  Now  rotate  these  new  axes  (that  is,  the  aircraft!)  138.5°  around  the  latest  snapshot 
of  ex  (north,  and  the  aircraft’s  nose),  which  brings  them  to  the  longitude  of  Adelaide, 
then 

3.  rotate  them  +34.9°  around  the  latest  snapshot  of  eY  (starboard  wing),  which  brings 
them  to  the  correct  latitude  of  Adelaide  and  now  makes  them  the  basis  vectors  of 
an  aircraft  flying  level  north  over  Adelaide,  then 

4.  rotate  the  axes  (the  aircraft)  45°  around  the  latest  snapshot  of  ez  (down),  which 
yaws  the  aircraft  to  head  north-east,  then  finally 

5.  rotate  the  axes  (the  aircraft)  20°  around  the  latest  snapshot  of  eY  (starboard  wing), 
which  pitches  the  aircraft  up  20°. 


This  sequence  of  rotations  is 


..Us  _  p 

t^E  —  A 


20° 


R 


45° 


(Y)  lh(Z)  1X-(Y) 


p+34.9°  pl38.5°  p— 90° 
-nW\  Hi  v\ 


L(X) 


1 Y 


/  T  \ 

[' ex\E  ieY]E  \ ez\E 


(5.24) 


Equation  (5.24)  is  a  recipe  for  constructing  the  basis  vectors  ex,ey,ez.  The  three  initial 
coordinate  vectors  [ex]E,  [ey\E,  [ ez]E  can  be  block  multiplied  by  the  relevant  rotation  ma¬ 
trices,  and  so  written  as  columns  of  a  matrix  in  (5.24) — which  is  just  the  identity  matrix, 
and  so  can  be  excluded  from  the  product: 

..Us  _  p20°  p45°  p+34.9°  pl38.5°  r-90°  /r  9^ 

HE  —  Hjyj  n(Z)  n(Y)  U(X)  UY  '  Z0J 

Now,  to  avoid  having  to  calculate  the  intermediate  “latest”  vectors  as  was  done  in  [1], 
we  invoke  the  Pseudo-Reversing  Theorem  to  convert  all  the  rotations  to  be  around  the 
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original  ECEF  axes,  thus  allowing  us  to  use  Euler  matrices  throughout: 


Ve 


o20°  o45° 

U(Y)  U{Z) 

p34.9°  o 
ix^y^  Jtx 

138.5° 

(X) 

O 

O) 

1 

fti 

0-90°  pl3S 
Ky  i;x 

.5°  p34.9° 

IXy 

Rf 

p20° 

lty 

0—9 0°  0138. 5°  o34.9° 
£/2 

Ef 

o 

o 

'-0.935 

-0.166 

0.313' 

-0.060 

-0.798 

-0.600 

0.349 

-0.580 

0.736 

(Pseudo- Reversing  theorem!) 


(5.26) 


The  required  orientation  matrix  is  the  transpose  of  this: 


-0.935 

-0.166 

0.313 


-0.060 

-0.798 

-0.600 


0.349 

-0.580 

0.736 


(5.27) 


Second  Way  to  Construct  :  This  approach  calculates  /i ^  directly.  It  doesn’t 
need  the  Pseudo- Reversing  Theorem,  and  is  the  same  sequence  described  in  the  First  Way 
above,  but  now  as  seen  by  the  aircraft  itself:  we  describe  events  from  our  viewpoint  in 
the  cockpit.  This  requires  us  to  imagine  that  Earth  (actually,  the  whole  universe)  rotates 
around  us.  We  now  rotate  snapshots  of  our  aircraft’s  x,  y,  z  axes  to  the  ECEF’s  X,  Y,  Z 
axes  by  doing  the  following. 

As  in  the  First  Way  above,  start  with  the  aircraft  initially  at  zero  lat/long,  with  its 
x,  y,  z  axes  parallel  to  the  ECEF’s  X,  Y,  Z  axes  respectively.  In  the  following  steps,  imagine 
rotating  Earth  to  bring  it  to  the  required  orientation  with  Adelaide  below  the  aircraft;  the 
aircraft  is  considered  to  be  at  rest  throughout  the  procedure.  Start  with  snapshots  of 
the  aircraft’s  basis  vectors  ex,ey,ez  which  will  be  embedded  in  Earth  itself  to  eventually 
become  the  required  ex,eY,ez  vectors.  Then  follow  the  motions  of  these  snapshots  as 
they  rotate  with  Earth  around  us. 


1.  Again,  the  aircraft  has  started  out  balanced  on  its  tail,  so  rotate  Earth  through  +90° 
around  the  aircraft’s  y  axis.  This  sets  the  ground  squarely  under  the  aircraft. 

2.  Now  rotate  Earth  by  —138.5°  around  the  aircraft’s  nose  ex,  which  brings  the  longi¬ 
tude  of  Adelaide  underneath  the  aircraft,  then 

3.  rotate  Earth  by  —34.9°  around  ey,  which  brings  Adelaide  exactly  under  the  aircraft, 
then 

4.  rotate  Earth  by  —45°  around  e2,  which  lines  north-east  up  ahead  through  the  wind¬ 
screen,  then 

5.  rotate  Earth  by  —20°  around  ey,  which  pitches  Earth’s  surface  down  by  20°. 


This  is  probably  harder  to  visualise  than  the  First  Way  above!  It  produces 


{ex,  •  •  • , ez}  —  R , 


-20°  ^>-45°  D-34.90  D-138.50  o90' 


y 


R 


y 


R , 


Ryu  ie  e  i 

±  *"U  I'-'X')  *  *  *  ?  '-'Z  /  ? 


(5.28) 
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so  that 


»U»  = 


[ ex\  Us  [er]  Us  i ez\  Us 


n-20°  n-450  0-34.9°  o-138.5°  o90° 
-tiy 


\ex\  Us  i^y\  Us  [®z]  Us 


-20°  f-45°  171—34.9°  ^-138.5°  J7i90' 


_  ZTi  —  ZU  iTi  —  40  77 

-  Uln  H/n  £/' 


1 


-^2  I 


=  identity  matrix 


(5.29) 


which  will  give  the  same  result  as  (5.26).  Why?  Because  to  invert  the  last  expression 
in  (5.29)  in  order  to  compare  it  with  (5.26),  we  must  reverse  the  order  of  its  matrix 
factors  and  invert  each  (that  is,  [ABCDE]^1  =  £'~1D”1C'_1I?_1A_1),  and  each  inverts 
by  changing  the  sign  of  its  rotation  angle. 

Whichever  way  we  choose  to  construct  /z|^,  we  can  now  write  (5.21)  as 


[rSA\us  -  V Us  irSA]E 


-0.935  -0.060 
-0.166  -0.798 
0.313  -0.600 


0.349 

-0.580 

0.736 


'-7.252' 

-9.213 

0.920 

x  105  metres 


7.654 

8.016 

3.933 


x  105  metres. 


(5.30) 


This  last  coordinate  vector  gives  the  other  aircraft’s  position  relative  to  us  in  our  coor¬ 
dinates:  the  other  aircraft  lies  at  a  point  roughly  765  km  ahead,  802  km  in  the  starboard 
direction,  and  393  km  down.  From  these  numbers  it’s  easy  to  see  that  we  sight  the  aircraft 
at  azimuth  (j>  (direction  positive  from  north  to  east)  and  elevation  6  (positive  upwards), 
where 

802  7  765  „  _i  -393  . 

sm  0  =  —  — 1  ,  cos(p= —  — -  .  v  =  tan  —  — —  .  5.31 

V7652  +  8022  V7652  +  8022  V7652  +  8022 

So  (j)  ~  46°  and  0  —  —20°:  we  in  the  cockpit  must  look  to  our  right  by  46°  and  down  20° 
to  sight  the  aircraft  flying  over  Sydney. 


5.3.1  A  Comment  on  Active  and  Passive  Rotations 

All  the  rotations  of  this  report  have  been  active ,  meaning  that  we  are  rotating  vectors.  The 
discussion  of  Section  5.3  above  is  essentially  worded  differently  in  some  books  on  rotation 
theory,  and  is  called  the  passive  rotation  approach.  The  idea  is  to  narrate  the  sequence 
of  rotations  of  the  First  Way  on  page  22  that  take  the  aircraft  from  its  datum  orientation 
at  zero  lat/long  and  produce  the  required  orientation  over  Adelaide,  while  accompanying 
this  with  a  merger  of  the  Second  Way’s  (5.29)  with  (5.30).  This  approach  allows  one  to 
ignore  the  fact  that  the  rotations  are  being  done  around  latest  axes  on  page  22,  in  that 
it  requires  a  sequence  of  rotations  to  always  use  latest  axes.  It  rewrites  each  of  these 
rotations  as  about  a  space- fixed  axis  through  minus  the  angle  turned  through,  giving  the 
Euler  matrices  in  (5.29) — which  has  the  effect  of  introducing  a  Ze/f-handedness  into  its 
conventions.  The  resulting  sequence  multiplies  (5.20)  to  give  the  required  [r\Us  in  (5.30). 
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Passive  rotations  are  all  about  rotating  a  set  of  axes,  in  contrast  to  the  active  approach’s 
explicitly  rotating  vectors.  Of  course,  a  set  of  axes  can  be  considered  as  a  set  of  basis 
vectors,  so  that  rotating  them  can  be  done  actively;  this  is  precisely  what  I  have  done  in 
this  report.  The  passive  approach  can  certainly  be  used  to  produce  (5.30).  But  I  think  it 
becomes  strained  when  we  really  wish  to  model  the  rotation  of  a  physical  object  described 
by  a  set  of  coordinate  vectors.  The  active  approach  that  I  have  used  throughout  this 
report  really  rotates  each  vector,  and  we  have  a  physical  picture  of  the  object  rotating.  In 
contrast,  the  passive  approach  doesn’t  explicitly  rotate  the  object.  Instead,  it  rotates  the 
coordinate  axes  the  other  way,  producing  all  coordinate  vectors  of  the  unrotated  object  in 
the  rotated  axis  set.  These  now  must  be  recognised  as  equalling  the  required  coordinate 
vectors  of  the  rotated  object  in  the  unrotated  axis  set.  Users  of  the  passive  approach  may 
or  may  not  be  aware  that  this  is  what  they  are  doing,  but  newcomers  to  the  subject  can 
find  it  all  very  confusing. 

The  active-rotation  approach  of  this  report  is  aimed  at  explicitly  constructing  coordi¬ 
nate  sets — of  which  there  might  be  several — and  then  forming  required  coordinate  vectors 
by  finding  components  using  the  dot  product.  This  is  all  neatly  encapsulated  in  (3.7). 


5.4  DIS  Conversion 

In  the  precursor  [1]  to  the  current  report,  I  gave  details  of  how  to  convert  position- 
orientation  information  into  and  out  of  the  IEEE  standard  [7]  known  as  the  Distributed 
Interactive  Simulation  environment,  or  DIS.  In  this  section  I  present  an  alternative  ap¬ 
proach  to  the  DIS  calculations  of  [1],  which  perhaps  will  better  disentangle  the  set  of 
calculations  necessary  to  implement  the  standard. 

DIS  position-orientation  information  of  an  object — say,  an  aircraft — is  recorded  in  the 
following  way.  The  aircraft’s  position  is  given  by  its  ECEF  ( X ,  Y,  Z )  coordinates.  Its 
orientation  is  specified  by  the  set  of  three  Euler  angles  that  enable  that  orientation  to  be 
constructed  by  three  rotations.  These  rotate  snapshots  of  the  ECEF’s  XYZ  axes  onto  the 
aircraft’s  xyz  axes.  The  three  rotations  are,  first,  by  the  angle  cp  around  the  ECEF’s  X, 
then  by  6  around  Y,  and  finally  by  ip  around  Z.  The  Pseudo-Reversing  Theorem  tells  us 
that  this  is  identical  to  rotating  first  by  ip  around  Z,  then  by  9  around  the  rotated  snapshot 
of  Y ,  and  finally  by  (p  around  the  latest  rotated  snapshot  of  X ;  this  is  the  sequence  as  it’s 
described  in  the  DIS  standard  itself,  so  that  the  three  angles  are  usually  specified  in  the 
order  of  ip,  9,  (p.  Writing  the  orientation  of  the  aircraft  A  relative  to  ECEF  E  as  /rg,  this 
all  means  that 

/4  =  EiE02Ef .  (5.32) 

Often,  the  orientation  of  the  aircraft  is  specified  by  its  location  along  with  three  non- 
DIS  Euler  angles:  the  angles  through  which  snapshots  of  the  local  NED  axes  would  be 
rotated  to  become  the  aircraft  axes.  How  do  we  convert  these  Euler  angles  to  the  DIS  Euler 
angles?  We  start  by  showing  how  to  extract  Euler  angles  from  any  orientation  matrix; 
to  be  completely  general  for  the  moment,  we  won’t  yet  insist  on  conforming  to  the  DIS 
convention  on  the  ranges  that  the  Euler  angles  must  take.  Remember  that  an  orientation 
matrix  is  a  rotation  matrix,  as  we  showed  in  (3.24).  In  fact,  it’s  not  hard  to  show  by 
construction  [4]  that  any  orientation  can  be  built  from  three  rotations  conforming  to  the 
DIS  standard,  although  we’ll  take  such  a  construction  as  a  given  here. 
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To  extract  the  DIS  Euler  angles  from  a  rotation  matrix,  let’s  write  a  general  rotation 
matrix  resulting  from  the  three  DIS  Euler-angle  rotations  above.  If  a  snapshot  of  axes  S 
is  rotated  through  cp,  9 ,  cp  in  the  DIS  way  to  produce  axes  S',  the  orientation  matrix  is 


^  =  E*  e‘2  E* 


cos  ip  —  sin  ip  0 
sin  cp  cos  ip  0 

0  0  1 


cos  # 

0 

sin# 

0 

1 

0 

sin# 

0 

cos  # 

1  0  0 

0  cos  cp  —  sin  (p 
0  sin  cp  cos  i p 


cos  ip  cos  6 
sin  ip  cos  9 
—  sin  9 


—  sin  0  cos  0 

-  cos  0  sin  6  sin  0 

cos  0  cos  0 

-  sin  0  sin  6  sin  0 


sin  0  sin  0 
+  cos  0  sin  6  cos  0 

—  cos  0  sin  0 
+  sin  0  sin  6  cos  0 


cos  9  sin  cp  cos  9  cos  (p 


(5.33) 


Call  this  matrix  R,  with  ijth  element  R^.  We  calculate  each  angle  by  stipulating  its  sine 
and  cosine. 


Remember  this  oft-neglected  point:  to  specify  an  angle  a,  we  always  need  two  pieces  of  trigono¬ 
metric  information.  In  general,  it’s  necessary  and  sufficient  to  specify  one  member  of  the  set 
{sin  a,  cos  a,  tan  a},  together  with  the  sign  of  another  member  of  this  set.  Usually  the  sine  and 
cosine  are  specified,  and  many  computer  languages  have  a  function  that  takes  these  two  numbers 
(or  any  positive  multiple  of  them)  as  arguments  and  returns  a.  This  function  is  typically  called 
atan2,  but  there  is  no  universally  agreed  order  for  its  arguments,  so  you  should  always  check.  Note 
that  atan2  is  a  computer  term  and  not  a  mathematical  function.  Many  people  make  the  mistake 
of  writing  a  =  tan-1  ,  but  the  function  tan-1  is  not  quite  defined  in  this  way,  and  it  will  only 
return  the  correct  angle  if  —90°  <  a  <  90°. 


If  cos  9  /  0  (i.e.  i?31  does  not  equal  1  or  —1),  then  we  can  see  by  inspecting  (5.33)  that 

•  ,  ^21  ,  -Rn 

smcp  = - -  ,  cos  ip  = 


cos  9  cos  9 

sin#  =  —  R3i  ,  cos#  =  sy/ 1  —  sin2  9 ' ,  where  s  £  {±1}; 


R 


sin  0  = 


■32 


R 


cos  9 


COS  0  = 


33 


cos  9 


(5.34) 


There  are  two  triplets  of  angles  here,  depending  on  an  arbitrary  choice  of  s.  Also,  if 
cos#  =  0  (i.e.  R3i  equals  1  or  —1),  then  sin#  can  be  1  or  —1  corresponding  to  #  =  ±90°: 


#  =  90°,  implying  R  = 

or  #  =  —90°,  implying  R  = 
In  other  words, 


0  —  sin(^>  —  cp)  cos  (ip  —  cp) 
0  cos  (ip  —  cp)  s m(ip  —  cp) 

-10  0 


0  —  sin  (ip  +  cp)  —  cos  (ip  +  cp) 
0  cos(cp  +  cp)  —  sin  (ip  +  cp) 

1  0  0 


(5.35) 


R3 1  —  +1  #  —  —90  ,  sin(t/i  +  cp)  —  —R12  —  —-^23:  cos(ip  +  cp)  —  — i?i3  —  R22  1 

R31  =  —1  ^=>  #  =  +90°,  sin  (ip  —  cp)  =  —R\2  =  +R231  cos  (V’  —  (p)  =  +^13  =  R22  ■  (5.36) 


Both  of  these  cases  yield  an  infinite  number  of  solutions  for  ip  and  cp,  but  they  are  all  valid; 
all  yield  the  same  orientation  of  the  aircraft. 
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For  DIS,  what  of  the  choice  of  sign  s  in  (5.34)?  The  two  choices  of  s  yield  two  sets  of 
Euler  angles,  and  if2i  $2 >02-  It’s  not  hard  to  show  that  each  set  is  as  good  as 

the  other,  although  we  won’t  do  that  here.  They  are  related  by 

V>2  =  180°  +  Vh  , 

e2  =  iso0  -  el , 

fa  =  180°  +  fa .  (5.37) 


The  DIS  standard  actually  requires 

-180°  ^  fa  (j)  ^  180°,  -90°  ^  0  <  90°.  (5.38) 

This  means  that  cos0  ^  0,  which  corresponds  to  choosing  s  =  1  in  (5.34).  If  i?31  does 
equal  1  or  —1,  we  can  solve  (5.36)  for  whatever  if  and  <f  obey  (5.38). 

This  last  case  of  0  =  ±90°  allows  for  multiple  choices  of  if  and  (f,  and  generally  any 
numerical  algorithm  that  computes  the  Euler  angles  will  become  difficult  to  manage  as  the 
aircraft  nears  that  orientation,  especially  if  the  algorithm  is  predicting  the  Euler  angles. 
Having  6  =  ±90°  corresponds  to  the  nose  pointing  toward  one  of  the  celestial  poles — 
generally  not  an  extreme  orientation  to  fly.'  For  example,  such  would  be  the  case  for  an 
aircraft  on  the  Equator  flying  level  and  pointing  north  or  south. 

The  fact  that  an  infinite  number  of  such  triplets  will  suffice  for  6  =  ±90°  means  that  Euler  angles 
can  present  difficulties  to  numerical  algorithms,  because  they  are  capable  of  changing  too  quickly 
for  any  algorithm  to  cope  with  when  9  nears  either  of  these  two  values.  As  a  result,  they  are  often 
viewed  with  suspicion  by  their  users,  as  if  something  mathematical  just  isn’t  quite  right.  But  Euler 
angles  do  exactly  what  they  have  been  set  up  to  do,  and  they  behave  completely  as  they  should 
once  we  understand  them,  even  though  this  behaviour  can  conflict  with  initial  expectations  that 
might  misinterpret  what  the  angles  were  designed  to  do. 

Note  that  our  extraction  of  the  Euler  angles  in  (5.34)  didn’t  use  all  of  the  matrix  /t|/.  This 
is  fine;  being  an  orientation  matrix,  (if!  has  internal  symmetry  that  means  its  elements 
are  not  all  independent  of  one  another.  The  omitted  elements  can  be  used  to  add  some 
numerical  strength  to  the  calculation,  but  such  a  study  in  numerical  stability  is  outside 
the  scope  of  this  report. 

With  the  extraction  of  the  DIS  Euler  angles  under  our  belt,  we  can  now  construct  DIS 
orientation  information.  We’ll  rework  the  example  of  [1]  here  differently  to  how  it  was 
done  in  that  report.  Here  is  the  original  question  posed  in  [1]: 

Convert  Lat  Long— Heading  Pitch— Roll  to  DIS  Position— Orientation  Coordi¬ 
nates:  An  aircraft  is  flying  10,000  m  above  Adelaide,  pointing  south-east,  climbing  at 
a  20° pitch,  and  holding  a  30°  roll.  What  are  the  two  sets  of  triplets  that  the  DIS  standard 
requires  to  specify  the  aircraft’s  location  and  orientation ? 

The  aircraft’s  DIS  location  is  its  ECEF  ( X ,  Y ,  Z)  coordinate  triplet.  This  is  found  by 
applying  (5.3)  to  Adelaide  (lat  —34.9°,  long  138.5°,  aircraft  height  10,000m),  resulting  in 
the  first  DIS  triplet 


( X,Y,Z )  =  (-3.928,  3.475,  -3.634)  x  106  nr.  (5.39) 

7I  thank  David  Sambcll  of  DSTO  for  pointing  this  out. 
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The  DIS  Euler  angles  ip,9,(j)  will  be  extracted  from  the  orientation  matrix  of  the 
aircraft  A  in  the  ECEF  E,  using  (5.32)  and  the  equations  immediately  after  it.  The 
ECEF’s  axes  are  X,Y,Z  and  the  aircraft’s  axes  are  x,y,z.  Remember  that  from  (3.11), 
for  example,  [ex]E  =  [ex\Ei  so  that  /i^  is  the  rotation  matrix  rotating  snapshots  of  the 
ECEF’s  basis  vectors  onto  the  aircraft’s  basis  vectors,  and  this  is  precisely  the  matrix  from 
which  we  must  extract  the  Euler  angles  using  (5.34);  these  will  be  the  DIS  Euler  angles. 

We  construct  from  the  following  sequence  of  rotations.  Align  the  aircraft  with  its 
x,  y,  z  axes  coincident  with  the  ECEF’s  X,  Y,  Z  axes.  Now: 

1.  rotate  snapshots  of  ex,eY,ez  around  eY  by  —90°, 

2.  rotate  the  resulting  three  vectors  around  the  rotated  version  of  ex  by  the  longitude 
u  =  138.5°, 

3.  rotate  the  resulting  three  vectors  around  the  rotated  version  of  eY  by  minus  the 
latitude,  —a  =  +34.9°, 

4.  rotate  the  resulting  three  vectors  around  the  rotated  version  of  ez  by  heading  135°, 

5.  rotate  the  resulting  three  vectors  around  the  rotated  version  of  eY  by  pitch  20°, 

6.  and  finally  rotate  the  resulting  three  vectors  around  the  rotated  version  of  ex  by 
roll  30°. 

This  is 


(Z)  ^(Y) 


'Y  J 


*4  (which  =  E$  E%  E*)  =  R^  fl]|f  B+if  >°  R^° 

which  the  Pseudo-Reversing  Theorem  converts  to 

He 


(5.40) 


A  {71—90°  {7il38.5°  t? 34.9°  eU35°  ci20°  ci30° 

—  jLVq  jLv  |  jL_/9  jC/q  iZ/9  JLj  \ 


J1 

-0.366  0.928  0.065 

-0.564  -0.165  -0.809 
-0.741  -0.333  0.584 


(5.41) 


The  (3, 1)  element  of  this  matrix  is  nonzero,  so  apply  (5.34) — with  s  =  +1  to  select  the 
DIS  convention — to  give  (retaining  more  decimal  places  internally  for  the  final  calculation 
of  the  angles) 


smi/j  = 


-0.564 


cos  9 

sin  9  =  0.741  , 

.  ,  -0.333 

sm0  = 


-0.366 

cos  = - x- ; 

cos  9 

cos  9  =  \/ 1  —  sin2  9 ' ; 
0.584 


cos  9 


COS  0  = 


cos  9 


(5.42) 


These  have  enough  information  to  yield  ( — refer  to  the  small  print  on  page  26) 

i/j  =  —122.97°,  9  =  47.79°,  (j)  =  -29.67°, 


(5.43) 
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which  are  the  required  DIS  angles.  Note  that  we  could  apply  (5.37)  to  produce  the 
alternative  set  of  Euler  angles  that  correspond  to  taking  s  =  —  1  in  (5.34): 

VLit  =  57.03°,  0alt  =  132.21°,  <^alt  =  150.33°.  (5.44) 

As  expected,  0alt  is  not  in  the  required  DIS  range  of  (5.38)  here.  (A  negative  cosine  will 
place  an  angle  into  the  range  90°  — ►  270°.  A  value  of  6  =  90°  is  allowed  by  DIS,  but 
cos  90°  =  0,  so  in  that  trivial  case  either  choice  of  s  can  be  used.)  So  the  angles  in  (5.44) 
don’t  conform  to  the  DIS  convention.  But  they  are  just  as  correct  as  those  of  (5.43)  as 
far  as  giving  the  aircraft’s  orientation  is  concerned. 

The  second  DIS  example  that  was  done  in  [1]  was  the  inverse  of  the  above  procedure, 
and  we’ll  repeat  it  in  greater  detail  now. 


Convert  the  Above  DIS  Position— Orientation  Coordinates  of  (5.39)  and  (5.43) 
Back  to  Lat— Long— Heading  Pitch— Roll:  An  aircraft  is  at  a  location  and  orientation 
given  by  the  above  DIS  coordinates 

(. X,Y,Z )  =  (-3.928,  3.475,  -3.634)  x  106  m, 

(VsM)  =  (-122.97°,  47.79°,  -29.67°).  (5.45) 


What  are  its  lat-long-height  on  Earth,  and  what  is  its  orientation  relative  to  local  north- 
east-down  in  terms  of  the  heading -pitch-roll  that  a  pilot  flies? 

(Again,  I  have  carried  through  more  decimal  places  from  above  and  have  rounded  them 
here.)  The  aircraft’s  lat-long-height  comes  from  (X,Y,Z)  and  was  calculated  in  [1],  but 
we  will  cover  some  of  the  ground  again  here.  We  require  to  solve  (5.3)  for  latitude  a , 
longitude  u,  and  height  h.  There  are  various  methods  for  solving  these  equations.  The 
approach  given  in  [1]  is  again  used  here:  first,  calculate  oj  from 


Y 


X 


sm  u>  = 


cos  u  = 


VX2  +  Y2  '  yj X2  +  Y2 

Now  iterate  to  find  a.  Begin  with  a  first  estimate 

'  f2  Z  \ 


(5.46) 


a  —  tan 


-i 


h 2v^ 


:2  +  y2 

with  a  and  b  from  (5.2).  Now  refine  this  estimate  using 


(5.47) 


a  =  tan 


-i 


9  •  2 

a  sm  a 


b2  sin  a  cos  a  +  ^  y/ X2  +  Y2  '  sin  a  —  Z  cos  y/ a 2  cos2  a  +  b2  sin2  a  f 


(5.48) 

(which  assumes  Z  0,  but  Z  =  0  is  a  trivial  case  anyway).  This  expression  for  a  converges 
very  quickly,  after  which  we  calculate  the  aircraft’s  height  h.  We  can  use  either  of 


h  = 


\! X2  +  y2 '  a2 

cos  a  \l a 2  cos2  a  +  b2  sin2  a ' 


(5.49) 
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or 

7  h2 

h  =  ~. - .  ,  (5-50) 

sm  a  cog2  q.  _|_  ft 2  sjn2  a 

depending  on  which  of  these  formulae  gives  the  best  numerical  stability;  that  is,  if  sin  a  ~  0, 
use  (5.49);  otherwise  use  (5.50).  Using  this  procedure  with  (X,Y,Z)  from  (5.45)  returns 
the  values 

latitude  a  =  —34.90°,  longitude  uj  =  138.50°,  height  h  =  10,000.00  m,  (5.51) 

meaning  10,000  metres  above  Adelaide,  as  expected. 

We  now  find  the  aircraft’s  heading- pitch-roll  relative  to  local  north-east-down.  These 
are  precisely  the  three  Euler  angles  that  would  be  required  to  rotate  snapshots  of  the 
north-east-down  axes  onto  the  aircraft’s  axes,  in  DIS  order.  They  will  be  extracted  from 
the  orientation  matrix  of  the  aircraft  relative  to  local  north-east-down  using  the  same 
approach  as  in  the  first  DIS  example  above.  The  orientation  matrix  is  where  A  is 
the  aircraft  and  N  is  the  NED  axis  set.  The  DIS  numbers  refer  to  the  ECEF  E,  so  first 
use  (3.9)  to  write 

/4  =  ^nVe=  ■  (5-52) 

We  calculate  /i^  from  the  just-found  latitude  and  longitude  (5.51),  and  calculate 
directly  from  the  DIS  ip,6,<p  (5.32).  We  showed  how  to  construct  /z^  in  Section  5.2,  but 
will  do  it  again  here  in  an  abbreviated  way.  Take  copies  of  the  ECEF’s  basis  vectors  and 
rotate  them  first  around  the  Y  axis  by  —90°,  then  around  the  latest  X  by  the  longitude, 
then  around  the  latest  Y  by  minus  the  latitude: 

,.N  p— a  TjuJ  p— 90° 

A lE  ~  n(Y)  U{X) 

=  E^90  E i  E2a  by  the  Pseudo-Reversing  Theorem,  (5.53) 

remembering  that  this  is  only  one  of  many  different  ways  in  which  we  could  construct  /z^. 
Equation  (5.32)  gives  /z^  =  E f  E®Ef,  so  (5.52)  becomes 

/4  =  (/4)_1a4 

=  E%  E^w  Ef°  E%  E°Ef 


-34.9°  rp- 
2  Eq 

138.5°  p9C 

E'2 

°  122.97°  p47.79°  p-29.67° 

-0.664 

— 0.733 

0.144' 

0.664 

— 0.491 

0.563 

-0.342 

0.470 

0.814 

The  heading,  pitch,  and  roll  are  extracted  from  this  matrix  by  our  realising  that  these  three 
angles  correspond,  respectively,  to  the  '0, 6,  that  would  construct  this  matrix  in  the  DIS 
rotation  order:  you  can  picture  the  process  as  the  pilot  beginning  with  nose  pointed  north, 
straight  and  level,  then  yawing  the  aircraft  onto  a  heading  of  ijj,  pitching  it  up  around  the 
starboard  wing  by  6,  then  rolling  it  about  its  nose  by  </>.  So  we  extract  9,  (f>  from  /.r(l  and 
relabel  them  heading,  pitch,  and  roll  respectively.  Be  careful  not  to  confuse  the  new  ip,  9,  cp 
here  with  the  ones  calculated  in  (5.42):  they  were  the  answer  to  a  different  question! 


30 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO-TR-2675 


The  (3, 1)  element  of  is  nonzero,  so  apply  (5.34)  with  s  =  +1: 


sin-0  = 


0.664 


cos  6 
sin  0  =  0.342 
.  ,  0.470 


sm  0  = 


-0.664 

cos-0  = - —  ; 

cos  t) 

cos  6  =  \J  1  —  sin2  O' ; 
,  0.814 

COS  0  = 


(5.55) 


cos  6  cos  6 

These  lead  to  the  expected  original  angles  (where,  as  usual,  I  have  used  more  decimal 
places  internally  and  then  rounded) 


heading  = -0  =  135.00°,  pitch  =  6  =  20.00°,  roll  =  0  =  30.00° 


(5.56) 


So  the  exercise  is  finished:  we  have  calculated  the  aircraft’s  position  in  (5.51)  and  its 
orientation  as  flown  by  a  pilot  in  (5.56).  We  also  have  its  orientation  matrix  in  (5.54), 
which  is  available  for  use  in  other  calculations. 


6  Dead-Reckoning  Aircraft  Orientation 


We  can  use  the  Pseudo-Reversing  Theorem  as  an  alternative  way  to  derive  the  rate  of 
change  of  an  aircraft’s  orientation,  given  its  angular  velocity  as  a  function  of  time.  This 
angular  velocity  is  typically  measured  by  gyroscopes  onboard  the  aircraft  that  make  all 
measurements  relative  to  its  body  axes.  The  aircraft’s  orientation  at  time  t  is  specified  by 
a  procedure  that  constructs  this  orientation  starting  from  some  base  orientation.  We  ask 
how  the  orientation  parameters  evolve  as  functions  of  the  angular  velocities  around  each 
of  the  body  axes,  as  measured  by  the  gyros.  We’ll  answer  the  question  for  three  different 
ways  of  quantifying  orientation:  the  orientation  matrix  in  Section  6.1,  angle-axis  using 
quaternions  in  Section  6.2,  and  Euler  angles  in  Section  6.3. 

For  this  analysis  it  proves  useful  to  show  that,  unlike  angular  displacement,  angular 
velocity  is  a  vector,  so  that  angular  velocities  can  be  added  in  any  order  to  give  a  combined 
angular  velocity.  This  can  be  proved  in  the  following  way.  First,  the  determination  of 
angular  velocity,  and  the  subsequent  integration  to  calculate  orientation,  is  a  first-order 
process,  meaning  that  exact  expressions  for  rates  of  increase  result  even  though  we  need 
“only”  calculate  to  order  d t.  This  first-order  approach  greatly  reduces  the  labour  of  the 
necessary  calculations,  but  why  does  it  work?  When  calculating  aircraft  orientation  over 
time,  we  can  work  with  time-dependent  quantities  by  using  Taylor’s  Theorem.  In  the 
more  familiar  language  of  position  s(t)  and  velocity  v(t)  in  one  dimension,  note  that 
velocity  is  a  first-order  increase  in  position  with  respect  to  time.  While  we  might  Taylor- 
expand  s(t  +  At)  to  write 


v(t)  =  lim 
v  ’  Ai— >o 


=  lim 

At^o 


s(t  +  At)  —  s(t ) 

At 

s(t)  +  s\t )  At  +  ^s"(t)  At2  + 


-  s{t) 


At 


1 


=  lim  s'(t)  +  — s"(t)  At  +  . . . 
At^o  w  2! 


=  s(t) 


(6.1) 
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we  can  write  this  analysis  in  precis  as 

=  »(t  +  dt)  -  »(t)  =  sjfl  +  s\t)  dt  -  »(t)  = 

w  dt  dt  w  v  ' 

Equation  (6.2)  is  a  shorthand  way  of  writing  (6.1).  Writing  “dt”  is  a  way  of  indicating 
that  we  are  effectively  writing  At  +  terms  of  order  At2  and  higher,  along  with  the  proviso 
that  at  some  stage  we  will  divide  by  At  and  take  a  limit  as  At  — >  0.  In  that  case,  terms 
of  order  At2  and  higher  won’t  survive  if  there  is  a  At  present,  and  so  we  needn’t  bother 
writing  those  higher-order  terms.  We  signal  that  we’re  following  this  procedure  by  writing 
the  At  as  dt  in  (6.2),  which  is  by  definition  exact. 

Position  is  found  by  integrating  velocity;  that  is,  multiplying  velocity  by  an  infinitesimal 
time  step  and  incrementing  the  previous  position  value, 

s(t)  =  s(tQ)  +  f  v(t)  dt,  (6.3) 

Jt0 

so  that  the  process  of  integrating  can  be  viewed  as  simple  multiplication  by  an  infinitesimal 
and  adding. 

For  infinitesimal  rotations  expressed  in  some  arbitrary  coordinates  S,  (3.14)  or  (3.28) 
both  give 

[i^]5  =  l  +  d@[n]x,  (6.4) 

where  “1”  is  the  unit  matrix,  and  we’ll  omit  the  [  ]g  notation  for  brevity  in  the  next  few 
lines.  As  noted  a  few  paragraphs  up,  the  way  that  infinitesimal  notation  is  defined  means 
that  this  is  an  exact  expression:  there  are  no  higher-order  terms  to  be  written.  The  key 
point  to  notice  is  that  we  are  now  able  to  represent  by  the  vector  d On,  because  this 
allows  a  double  rotation  R ^  R ^  to  be  written  as8 

R*%  Rf  =  (1  +  damx)(l  +  d/3nx) 

=  1  +  da  mx  +  d/3  nx 

=  1  +  (da  m  +  d(3  n) x ,  (6.6) 


which  can  thus  be  represented  by  da  m  +  d/3  n,  as  required  if  the  vector  representation  is 
to  be  meaningful.  So  infinitesimal  rotations  behave  as  vectors,  and  combining  rotations  is 
represented  by  adding  their  vectors.  Non-infinitesimal  rotations  don’t  behave  as  vectors 
because  they  don’t  commute,  so  they  don’t  conform  to  the  behaviour  of  vector  addition, 
which  is  commutative.  This  corresponds  to  attempting  to  use  A  On  plus  higher-order 
terms  in  (6.6):  those  higher-order  terms  will  prevent  any  identification  of  the  rotation 
with  A  On. 

We’re  now  in  a  position  to  add  angular  velocities  by  examing  the  effect  of  each  over  an 
infinitesimal  time  dt.  Dividing  the  corresponding  infinitesimal  rotations  d@x  rii  and  d 02  n2 
by  dt  retains  their  vector  character;  but  this  just  produces  vectors  whose  lengths  are  the 
corresponding  angular  velocities: 


df?!  d  02 

n1;  u2  =  — n2 


(6.7) 


8It’s  useful  in  (6.6)  to  remember  that  the  operation  (3.15)  of  “applying  the  cross”  to  a  vector  is  linear, 
meaning  that  for  scalars  a,  /3, 


(am  +  f3n)x  =  amx  +  j3nx . 


(6.5) 
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So  angular  velocity  is  indeed  a  vector.  Given  tiq  and  u>2,  over  a  time  d t  their  combined 
infinitesimal  rotation  is  represented  by 


d$!  n1  +  d d2  n2  =  uq  dt  +  uj2  d t 
=  (uq  +  lj2)  dt , 


(6.8) 

(6.9) 


so  that  the  combined  angular  velocity  is  +  cj2,  as  expected. 

With  this  in  mind,  consider  the  three  angular  velocities  measured  by  the  gyros  onboard 
an  aircraft.  These  make  their  measurements  in  an  inertial  frame  that  momentarily  shares 
the  aircraft’s  linear  velocity.  They  report  that  the  aircraft’s  angular  velocity  about  its 
x  axis  (nose)  is  p(t)  (we’ll  omit  the  t  for  brevity),  about  its  y  axis  (starboard  wing) 
is  q(t),  and  about  its  z  axis  (the  below-fuselage  direction)  is  r(f).  In  body  coordinates  B, 
the  angular  velocity  coordinate  vector  about  the  nose  is  thus  (p,  0,  0),  the  angular  velocity 
coordinate  vector  about  the  starboard  wing  is  (0,  q,  0),  and  the  angular  velocity  coordinate 
vector  about  the  below-fuselage  direction  is  (0,0,  r).  The  combined  angular  velocity  must 
then  be  the  sum  of  these: 


Mb  - 


P 

q 


(6.10) 


r 


The  Pseudo-Reversing  Theorem  clears  up  any  reservations  we  might  have  about  how  the 
angular  velocities  are  to  be  understood,  for  suppose  that  these  velocities  are  visualised 
as  rotations  that  act  first  through  angle  pdt  around  the  nose  ex,  then  through  angle  qdt 
around  the  latest  wing  e^\,  then  through  angle  rdt  around  the  latest  below-fuselage 
direction  e^zy. 

Rpxdt  ^  Rr$ .  (6.11) 

The  PRT  says  that  this  is  equivalent  to  R^dt  Rydt  Rrzdt.  We  see  that  any  sense  of  “latest” 
axes  has  now  been  discarded,  and  furthermore,  these  three  infinitesimal  rotations  com¬ 
mute,  so  can  be  applied  in  any  order.  This  is  just  as  well,  since  the  gyros  do  make  their 
measurements  simultaneously  and  so  no  rotation  is  necessarily  “later”  than  any  other. 

We’re  now  in  a  position  to  calculate  the  rate  of  change  of  an  aircraft’s  orientation  with 
time,  as  a  function  of  p(t),q(t),r(t). 


6.1  Calculating  Rate  of  Change  of  an  Orientation  Matrix 

The  aircraft’s  orientation  can  be  specified  by  the  matrix  p^(t)  whose  columns  are  the 
aircraft  body  ( B )  basis  vectors  specified  in  some  world  coordinates  W.  Write  the  body 
axes  as  x,  y,  z.  The  body  basis  vectors  are  rotated  over  a  time  interval  dt  by  the  angular 
velocity  u  =  urn,  so  that 


—  \ex{t  +  dt)]jy  .  .  .  [e2(t  +  dt)]jy 


—  [Rndt]w  [ex(t)]w  ■  ■  ■  [ez{t)\w  —  (l  +  [w]w  dtj  ■  (6.12) 


We  conclude  that  the  time  derivative  of  the  orientation  matrix  is 

,B  u  i  ,.B 


/£(*)  = 


p,w(t  +  dt)  -  pw(t) 
dt 


=  u? 


PwW- 


w 


(6.13) 
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But,  referring  to  (6.10),  the  aircraft  gyroscopes  give  us  [u>]#,  not  these  two  coordi¬ 

nate  vectors  are  related  by  the  orientation  matrix  /i^,,  but  this  matrix  is  precisely  what 
we  are  trying  to  calculate.  Address  this  by  interchanging  B  and  W  in  (6.13),  remember¬ 
ing  that  reversing  the  viewpoint  requires  u>  — >  —  u>,  since  the  aircraft  sees  the  world  spin 
oppositely  to  how  the  world  sees  the  aircraft  spin: 

=  (6.14) 

But  we  require  (i^v ,  not  /i^  !  Remedy  this  by  noting  that  p}'^  p^v  =  1,  which  differentiates 
to  =  0,  or 

Pb  V'W  ~  ^ B  P\v  ’  (6.15) 

That  is,  we  can  shift  the  time  derivative  from  one  factor  to  the  other  in  (6.15)  at  the  cost 
of  a  minus  sign.  Now  post-multiply  (6.14)  by  to  give 


(6.16) 

which  (6.15)  changes  to 

(6.17) 

Premultiplying  this  by  gives 

Pw  =  /4  M  B  ■ 

(6.18) 

[Compare  this  with  (6.13).]  Equation  (6.18)  is  the  standard  expression  for  the  time  evolu¬ 
tion  of  the  aircraft’s  orientation  matrix,  using  the  gyro-supplied  rotation  rates  p,  q ,  r.  That 
is,  an  inertial  navigation  system  can  take  these  rotation  rates  and  use  them  to  update  its 
knowledge  of  the  aircraft’s  orientation  over  time.  An  alternative  derivation  of  (6.18)  can  be 
found  in  [8].  (Remember  that  the  orientation  matrix  is  usually  called  the  direction-cosine 
matrix  in  the  literature.)  Relatively  long  time  intervals  between  gyro  updates  can  be  ac¬ 
commodated  by  solving  (6.18)  using  an  efficient  differential-equation  solver  that  maintains 
long-term  accuracy. 

6.1.1  An  Alternative  Derivation  of  Equation  (6.18) 

Equation  (6.18)  can  be  produced  in  a  more  generic  way  with  a  careful  application  of  the 
Pseudo- Reversing  Theorem.  Suppose  A(t)  is  some  mathematical  object  that  quantifies  the 
orientation  of  the  aircraft  at  time  f  by  specifying  how  to  rotate  the  world  basis  vectors  W 
to  the  body  basis  vectors  B.  In  practice  the  world  basis  vectors  are  all  we  have  to  work 
with,  so  we  require  [A(f)]^.  It  might  be  the  orientation  matrix  p^n  but  it  could  also  be 
something  more  exotic  and  novel.  We  require  the  dependence  of  [A(t)\w  on  p,  q,r. 

By  definition,  the  aircraft’s  orientation  at  time  f  is  produced  by  having  A(t)  act  on  a 
base  orientation.  In  the  next  time  step  df,  the  aircraft  rotates  around  u  =  pex  +  qey  +  rez 
(remember  that  ex,ey,ez  are  body  basis  vectors),  and  the  new  orientation  is  specified 
by  A(t  +  df): 

A(t  +  df)  =  A(t)  — >  rotate  around  u  =  pex  +  qey  +  re~  by  u  df  .  (6.19) 
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The  PRT  converts  this  to  rotations  around  the  W  axes  xw,yw,  zw: 

A(t  +  d t)  =  A(t)  (rotate  around  peXw  +  Qeyw  +  rezw  by  u>  dt)  .  (6.20) 

In  W  coordinates  this  is 

[A(t  +  dt)]w  =  [A(t)\w  [rotate  around  peXw  +  qeyw  +  reZw  by  udt]w.  (6.21) 

The  theorem  doesn’t  alter  the  sense  of  rotation  of  A(t),  as  you  can  see  by  recalling  (2.10). 
Now,  the  key  point  to  notice  is  that 

[Pexw  +  qeyw  +  reZw ]  w  =  [pex  +  qey  +  rez\  B  (6.22) 

r  p~\ 

since  both  equal  q  .  This  converts  (6.21)  to 

r 

[A(t  +  df)]^/  =  [A{t)]w  [rotate  around  u>  by  udt]B  .  (6.23) 

For  example,  if  [A(t)\w  is  chosen  to  be  then  (6.23)  becomes 

P-wtt  +  dt)  =  (x  +  Ms  dt)  ■  (6-24) 


This  leads  very  quickly  to  (6.18),  as  expected.  Equation  (6.23)  is  a  generic  equation  that 
can  be  applied  to  any  choice  of  [A(t)]w,  as  we’ll  see  next  when  we  use  a  quaternion  to 
quantify  the  aircraft’s  orientation. 


6.2  Calculating  Rate  of  Change  of  Quaternions  for  Angle- 
Axis  Representation 

Set  [A(t)\w  to  be  the  quaternion  Q^(t)  that  rotates  the  W  basis  to  the  B  basis.  Equa¬ 
tion  (6.23)  becomes  a  quaternion  multiplication 

Qw(t  +  di)  =  Qwtt)  Q[n]fB  (6.25) 

(where  m  =  can),  because  combining  multiple  rotations  using  quaternions  is  accomplished 
by  multiplying  the  quaternions,  as  we  saw  in  (3.35).  Now  refer  to  (3.30)  to  write  the 
quaternion  for  the  infinitesimal  rotation  as 

Qf4fl  =  (hMfidt/2)  =  (l,\p,q,r]dt/2)  =  (1,0, 0,0)  +  (0,p,  q,  r)  dt/2  .  (6.26) 

identity  quaternion  =  Q(t) 

Defining  a  body  angular  velocity  quaternion  f l(t)  =  (0 ,p,  q,  r)  (that  uses  the  angular  veloci¬ 
ties  in  body  coordinates)  allows  (6.25)  to  be  written  using  the  identity  quaternion  (1, 0,  0, 0) 
(the  quaternion  that  gives  zero  rotation,  which  can  just  be  written  as  “1”): 

Qw(t  dt)  =  Qw(t)  [l  +  n(t)  dt/2] 

=  Qw(t)  Qw(t)  ^(t)  dt/2  .  (6.27) 


This  reduces  to 

Qw(t)  =  2  Qw(t)  ^(0  • 

An  alternative  derivation  of  this  well-known  differential  equation  is  in  [8]. 


(6.28) 
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6.3  Calculating  Rate  of  Change  of  Euler  Angles  using 
Matrices 


The  aircraft’s  orientation  can  be  specified  by  giving  a  sequence  of  three  rotations  through 
Euler  angles  that  rotate  the  base  orientation.  These  rotations  are  conventionally  taken 
to  be  around  either  two  or  three  different  cartesian  axes,  [^(f)]^  is  usually  expressed  as 
a  matrix  when  Euler  angles  are  used;  sometimes  it’s  written  as  a  quaternion.  We’ll  use 
matrix  language  in  what  follows. 

Choosing  x,  y,  z  to  refer  to  e.g.  ECEF  axes  or  some  other  base  set,  and  invoking  the 
DIS  convention  on  page  25,  the  rotation  order  becomes 

A(t)  =  Rt(t)  ->  Re$  -  R$.  (6.29) 

(We’ll  omit  the  explicit  t  dependence  of  ip,  9 , cp.)  The  PRT  converts  this  to 

A(t)  =  RfReyRt,  (6.30) 

so  that 

[A(t)]w  =  EtEe2Ef, 

where  we  have  replaced  the  generic  rotations  with  the  Euler 
tion  (6.23)  now  becomes 

E|+dV’  Ee2+Ae  Ef+d<l>  =  Ef  Ee2  Ef  (l  +  [w]  *  dt) .  (6.32) 


(6.31) 

matrices  of  (3.19).  Equa- 


This  equation  must  be  solved  for  the  time  derivatives  ip,9,<p  (where  e.g.  ip  =  dip/dt)  in 
terms  of  p,q,r.  A  straightforward  but  tedious  approach  inserts  the  Euler  matrices  (3.19) 
into  (6.32)  and  then  multiplies  the  matrices  to  extract  ip,  9,  <p.  Much  easier  is  to  use  (6.4) 
to  write 


E^+di,  =  ^ 0  Ep>  =  E 0 


1  +  d  ip 


(6.33) 


and  similarly  for  the  other  matrices  of  the  left-hand  side  of  (6.32).  These  expressions 
allow  £3’  E2  Ef  to  cancel  from  both  sides  of  (6.32),  greatly  reducing  further  labour.  The 
remaining  manipulations  are  straightforward  and  yield  (as  in  [8],  which  uses  an  alternative 
approach) 


sec  9  sin  cp 

V 

U 

sec  9  cos  (p 

P 

9 

= 

0 

cos  cp 

—  sin  <p 

Q 

(p 

1 

tan  9  sin  <p 

tan  9  cos  <p 

r 

(6.34) 


This  way  of  calculating  ip,  9,  cp  easily  accommodates  a  change  to  the  Euler  order.  Suppose 
that  convention  (6.29)  was  changed  to  a  rotation  about  two  space- fixed  axes  with  one 
repeated,  defining  a  different  set  of  angles  a,  (3, 7: 


A(t)  =  R*  -»  Rf  RJ  .  (6.35) 

The  required  matrix  multiplication  is  immediately  written 

[A(t)]w  =  E%E$E%  (6.36) 
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without  needing  the  PRT,  after  which  the  calculations  proceed  in  the  same  way  as  those 
following  (6.31) — although  the  final  result  need  not  resemble  (6.34). 

Newcomers  to  the  subject  of  aerospace  orientation  sometimes  confuse  the  angular  ve¬ 
locities  p,  q,  r  measured  by  gyroscopes  with  the  Euler-angle  rates  of  increase  <p ,  6 ,  ip.  They 
ask  why  the  matrix  in  (6.34)  is  so  complicated;  shouldn’t  d <p  just  equal  pdt,  so  that  cp  =  p, 
and  similarly  for  0,ip7 

Repeating  some  of  the  previous  analysis  in  different  language  makes  the  situation  clear. 
Using  the  DIS  Euler  order,  the  orientation  at  time  t  is  produced  by  having  E f  E 2  Ef  act 
on  the  base  orientation.  In  the  next  time  step  dt,  the  aircraft  simultaneously  rolls  about 
its  nose  (the  latest  x)  by  angle  pdt,  pitches  about  the  latest  y  by  qdt,  and  yaws  about  the 
latest  z  by  rdt.  Using  the  PRT,  the  new  orientation — equation  (6.32) — becomes 

Ef  E^  Ee2  E2e  EfE^  =  Ef  E62  Ef  Er3dt  Eq2dt  E[dt .  (6.37) 

Both  sides  of  this  equation  contain  the  Euler  rotations  Ef ,  E2,  Ef;  but  whereas  its  left- 
hand  side  interleaves  these  rotations  with  infinitesimal  rotations  such  as  Ed^ ,  its  right- 
hand  side  implements  all  of  its  infinitesimal  rotations  Ef dt ,  E2  dt ,  Er3  dt  first,  before  doing 
the  three  Euler  rotations.  But  infinitesimal  rotations  commute  only  with  each  other  and 
not  with  non-infinitesimal  ones  (see  page  10),  so  we  are  obliged  to  retain  the  given  order 
throughout  (6.37).  This  implies  that  the  infinitesimal  angles  are  generally  different  on  each 
side  of  that  equation,  so  that  d <p  does  not  simply  equal  pdt.  The  same  idea  applies  to  dd 
with  qdt,  and  dip  with  rdt.  And  that  explains  why  the  matrix  in  (6.34)  is  so  complicated. 


7  Concluding  Comments 

The  use  of  the  Pseudo-Reversing  Theorem  is  part  of  a  larger  context  within  which  vectors 
exist  independently  of  coordinate  systems,  and  yet  must  be  expressed  using  coordinates 
for  use  in  numerical  calculations.  Such  expressions  require  the  machinery  of  coordinate 
changes.  We  are  free  to  calculate  using  any  coordinates  we  choose,  but  rotations  that 
are  easy  to  describe  and  implement  by  an  aircraft  pilot  are  not  always  easy  to  implement 
mathematically,  and  vice  versa.  The  PRT  connects  these  two  rotation  schemes.  It  validates 
choices  of  coordinates  and  rotation  order  that  can  sometimes  appear  mysterious  in  the 
literature  when  they  are  accompanied  by  little  or  no  explanation. 
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Appendix  A  An  Application  of  Orientation 
Matrices  to  Measured  Wind  in  Ship  Motion 

Section  3  of  this  report  stressed  the  distinction  between  proper  vectors  and  coordinate 
vectors,  and  used  (3.7)  to  convert  one  coordinate  representation  to  another.  This  view 
of  vectors  is  very  useful  for  untangling  the  concepts  involved  with  the  wind  that  sailors 
feel  to  be  flowing  over  their  ship.  In  this  appendix  we’ll  use  that  notation  to  answer  the 
following  often-asked  question: 

Calculating  a  Desired  Ship  Course:  A  ship  is  sailing  at  some  given  speed  on  some 
given  compass  bearing.  It  measures  the  speed  and  direction  of  the  wind  flowing  over  it. 
What  new  course  should  it  sail  so  that  the  new  speed  and  direction  of  the  wind  over  the 
ship  become  some  given  values? 

Use  the  following  notation. 

The  compass  bearing  is  always  relative  to  local  north-east-down,  so  let  N  stand  for 
these  coordinates  and  for  the  general  reference  frame  attached  to  these:  e.g.,  think 
of  N  as  also  representing  some  nearby  island. 

-  W  stands  for  the  wind. 

-  S  stands  for  the  current  ship  and  its  coordinates.  We’ll  take  the  x  axis  to  point 
forward  of  the  prow,  and  the  y  axis  to  point  to  starboard. 

-  S'  stands  for  the  “new  ship”  and  its  coordinates:  those  of  the  new  orientation  that 
we  require  to  sail.  The  axes  of  this  “new  ship”  will  be  just  as  for  the  current  ship, 
but  labelled  with  primes:  x'  pointing  forward  of  the  prow  and  y'  pointing  starboard. 

-  vAB  is  the  velocity  of  object  A  relative  to  object  B,  and  is  a  proper  vector.  So  for 
objects  A ,  B,  C, 

VAB  =  VAC  +  VCB  ■  (Al) 


It’s  worthwhile  covering  some  established  terminology  relating  to  the  wind. 

True  wind  is  vWN,  the  actual  motion  of  the  wind  over  Earth’s  surface,  and  is  a  proper 
vector.  We  can  think  of  it  as  the  motion  of  the  wind  relative  to  anything  at  rest  in 
the  local  north-east-down  system,  such  as  a  nearby  island. 

Measured  wind  is  vWB,  the  velocity  of  the  wind  relative  to  the  ship:  a  proper  vector. 

Apparent  wind  is  [i>u/5,]a''>  the  measured  wind  in  north-east-down  coordinates:  a  coor¬ 
dinate  vector. 

Relative  wind  is  [i>ws]s>  the  measured  wind  in  current  ship  coordinates:  a  coordinate 
vector. 

Note  that  there  are  only  two  proper  wind  vectors  here:  the  true  wind  (velocity  relative 

to  an  island)  and  the  measured  wind  (velocity  relative  to  the  ship).  The  apparent  and 
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relative  winds  are  coordinate  vectors:  just  two  different  ways  of  quantifying  the  measured 
wind. 

Now,  both  the  “current  ship”  and  “new  ship”  agree  on  what  the  true  wind  is  (assuming 
it  didn’t  change  during  the  manoeuvre,  of  course),  because  the  true  wind  is  a  proper  vector 
and  so  has  a  reality  independent  of  any  ship.  So  we  can  write  the  true  wind  vector  as  a 
sum  of  terms  using  the  current  ship  or  the  new  ship,  using  (Al): 

ywNt  =  ywsi  +  ysNi  =  yWs\  +  ys,Ni  ■  (A2) 

true  wind  current  wind  current  ship  desired  wind  new  ship 
over  ship  velocity  over  ship  velocity 

To  answer  the  main  question  posed  above,  we  seek  the  new  ship  course  vS'N,  which  comes 
from  (A2): 

VS'N  =  vws  +  VSN  ~  vws’  ■  (A3) 

This  is  the  basic  vector  sum  that  will  be  coordinatised  with  respect  to  either  north-east- 
down  or  the  ship,  depending  on  what  information  we  have  been  given.  The  main  question 
above  requires  us  to  coordinatise  it  as  follows,  using  (3.7): 


Its'jv ] N  -  [ vws]n  +  [vsn]n  -  [ vWS']n 

=  Mjv  JTwsIs  +  few]iV  -  ^7V  J^Wsfe  •  (A4) 

given  given  given 

The  first  “given”  term  in  (A4)  is  the  current  relative  wind,  which  is  supplied  in  the  question. 
The  second  “given”  term  is  the  current  ship  course.  The  third  “given”  term  is  the  desired 
relative  wind. 

Equation  (A4)  is  not  straightforward  to  solve,  because  the  desired  ship  course  on  the 
left-hand  side  of  (A4)  establishes  the  new  orientation  via  S',  but  this  orientation  is  required 
in  the  last  term  on  the  right-hand  side  of  (A4).  However,  we  can  find  an  exact  solution. 

Begin  by  specifying  all  coordinates.  We  who  stand  on  the  ship  measure  the  wind 
coming  at  us  with  speed  tq  from  an  azimuth  Pi  measured  clockwise  from  the  prow: 


Nds  = 


— Vi  COS  Pi 
—Vi  sin  /Sj 


(A5) 


The  velocity  of  our 
bearing  p2: 


ship  in  the  north-east-down  coordinates  is  speed  v2  along  compass 


few]  AT 


V2  COS  P2 
v2  sin/32 


(A6) 


We  require  to  turn  such  that  we’ll  measure  the  wind  to  be  coming  at  us  with  speed  v3 
from  an  azimuth  /?3  measured  clockwise  from  the  prow: 


[vws']s' 


-v3  cos/33 
-v3  sin/33 


(A7) 


The  sought-after  new  ship  velocity  is  described  by  some  number  v  along  compass  bearing  /?: 


fe'wlw  — 


v  cos  P 
v  sin  P 


(AS) 
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There  is  nothing  to  prevent  v  from  being  negative,  so  we’ll  refer  to  it  as  a  velocity  rather 
than  as  a  speed.  For  shorthand,  set  for  i  =  1,2,  3, 


c  =  cos  f3  ,  s  =  sin  (3  , 

q  =  cos  /3j ,  s,  =  sin  f3i .  (A9) 

Equation  (A4)  requires  the  orientation  matrices  /i)(,  and  .  We  need  to  know  the  ship’s 
orientation,  so  will  assume  it  moves  without  side-slipping.  The  notation  can  handle  any 
orientation  at  all,  but  ships  do  usually  move  this  way.  If  the  ship  currently  points  into  the 
bearing  direction  /32,  then 


/4 


c2  ~s2 
s2  c2 


Likewise,  we’ll  assume  the  “new  ship”  points  into  the  bearing  direction  (3: 


/4 


e„/ 


x'iN 


e„/ 


y'\N 


C  —s 
s  c 


With  these,  (A4)  becomes 


v 


L 


~s2 

C2 


Cl 

«1 


+  V2 


C2 

s2 


+ 


call  this 

a 

=  v3 

C3 

~s3 

C 

b 

s3 

c3  _ 

S 

which  can  be  rewritten  as 


c3 

«3 


~s3 

C3 


c 

s 


(A10) 


(All) 


(A12) 


(A13) 


where  by  “v.l”  we  mean  v  times  the  2x2  identity  matrix.  This  last  equation  becomes 


v  -  V3C3  V3S3 

c 

a 

-v3s3  v-v3c3_ 

s 

b 

(A14) 


This  must  be  solved  for  v,  c,  and  s.  To  do  so,  equate  the  squared  lengths  of  each 
side  of  (A14),  remembering  that  for  a  matrix  A  and  a  column  (coordinate  vector)  x, 
|Ax|2  =  xtAtAx.  (Remember  too  that  the  elements  of  AfA  are  just  dot  products  of  the 
columns  of  A  with  themselves  in  turn.)  Writing  the  squared  length  of  [&]  as  k 2,  we  obtain 


r  n  [  (v  -  v3c3)2  +  vl  s|  0 

C  s  n9  9  9 

L  0  (v  -  V3cs)  +  V3  4 

In  other  words, 

(v  —  v3c3)2  +  v3  s3  =  k2 ,  (A16) 

which  is  easily  rewritten  to  give  the  desired  new  ship  velocity  as 

v  =  v3c3±  sjk2-v l  a§'.  (A17) 


=  k2.  (A15) 
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There  could  be  any  of  zero,  one  or  two  possible  real  values  of  v  resulting  here.  Now  we 
need  only  use  whatever  values  of  v  we  have  obtained  here  to  invert  (A14),  making  use 
of  (A16): 


c 

v  -  U3C3 

Ui®3 

-1 

a 

s 

-v3s3 

v  -  U3C3  _ 

b 

1 

¥ 


V  -  U3C3 

US' S3 


-Us  Us 

V  -  U3C3 


(A18) 


This  last  equation  gives  us  cos  / 3  and  sin  0  and  hence  (3  itself,  remembering  the  comment 
on  page  26  that  two  pieces  of  trigonometric  information  are  always  needed  to  specify  an 
angle.  So  the  calculation  is  finished. 

As  an  example,  suppose  we  on  the  ship  measure  the  wind  coming  at  us  with  speed 
15  units  (it  doesn’t  matter  which  units  we  use,  as  long  as  they’re  all  the  same)  from  an 
azimuth  140°  measured  clockwise  from  the  prow: 


=  15,  ft  =  140°. 


(A19) 


Our  ship  has  speed  20  over  the  ocean,  along  compass  bearing  270°: 


u2  =  20,  02  =  27Oo. 


(A20) 


We  require  to  turn  such  that  we’ll  measure  the  wind  to  be  coming  at  us  with  speed  10 
from  an  azimuth  44°  measured  clockwise  from  the  new  position  of  the  prow: 


*3  =  10,  03=44°.  (A21) 

The  required  ship  velocity  is  v  along  compass  bearing  0.  Use  (A12)  to  calculate  [£],  then 
find  its  squared  length  k2,  then  apply  (A17)  to  produce  v.  Finally  calculate  0  from  its 
cosine  and  sine  in  (A18).  Two  pairs  of  u,0  result: 


v  =  39.4 ,  0  =  265.2°, 

(A22) 

v  =  —25.0  ,  0  =  60.8°. 

(A23) 

The  first  pair  describes  the  ship  moving  with  speed  39.4  into  compass  direction  265.2°, 
and  is  one  solution  to  the  original  question  posed. 

The  second  pair  cannot  be  read  as  describing  a  speed  of  25.0  into  a  compass  direction 
180°  +  60.8°  =  240.8°.  Rather,  we  used  the  new  direction  of  ship  travel  to  specify  its 
orientation  in  (All):  we  assumed  that  the  ship  points  in  the  direction  of  its  velocity,  or 
compass  direction  0.  That  means  the  second  pair  of  numbers  (A23)  must  describe  a  ship 
travelling  backwards  with  speed  25,  while  pointing  in  the  compass  direction  60.8°.  Such 
a  ship  would  indeed  feel  the  wind  coming  with  speed  10  from  an  azimuth  44°  measured 
clockwise  from  its  prow;  but  no  ship  is  going  to  move  in  reverse  like  this.  Even  so,  (A23) 
is  a  valid  solution  to  the  question  posed. 
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